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the  daily  load  of  events  which  must  be  subjected  to  classification  processing. 
Flexibility  in  altering  operating  procedures  will  be  needed  until  our  discrimi- 
nation capability  is  sufficiently  developed  to  identify  explosions  originating 
from  any  region  which  is  possibly  accessible  to  nuclear  explosions. 

To  partially  fulfill  these  requirements,  a special  purpose  Interactive 
Seismic  Processing  Module  was  developed  on  a PDP-15  minicomputer  for 
Short-Period  Earthquake/Fxplosion  Discrimination  (SPEED).  This  processing 
module  was  imbedded  in  a simulated  special  purpose  seismic  operating  sys- 
tem, Interactive  Seismic  Processing  System  (ISPS)  developed  by  Ringdal  and 
Shaub.  All  of  the  capabilities  for  waveform  processing  using  SPEED  were 
described  and  demonstrated  by  processing  a seismic  event  and  displaying  the 
results  of  each  stage  of  processing.  Support  programs  were  developed  to  use 
an  IBM  360/44  computer  to  reduce  all  of  the  station  component  seismograms 
to  a single  representative  event  waveform.  An  option  is  provided  for  input- 
ting into  SPEED  the  single  representative  event  waveform  either  corrected 
for  or  uncorrected  for  system  response.  The  examples  analyzed  were 
corrected  to  broadband  ground  displacement  between  0.  25  and  5.  0 Hz. 

The  SPEED  processing  module  as  presently  configured  can  be  pro- 
grammed to  perform  any  one  of  many  discrimination  procedures  by  a pre- 
scribed sequence  of  button  pushing.  To  obtain  baseline  performance  data, 
one  very  simple  and  reproducible  discrimination  processing  procedure  was 
prescribed.  A data  base  of  35  events,  including  20  earthquakes  and  15  pre- 
sumed explosions  in  the  magnitude  range  between  4.  4 and  6.  1,  were  run 
through  the  SPEED  processor.  For  each  event,  the  discrimination  measure- 
ments obtained  from  SPEED  were  processed  by  multivariate  discrimination 
analysis  of  frequency  dependent  magnitude  measurements  to  derive  a two- 
component  discriminant  output.  One  component  was  optimized  to  discrimi- 
nate central  Asian  presumed  explosions.  The  other  component  w’as  optimized 
to  discriminate  western  United  States  presumed  explosions.  The  perform- 
ance results  obtained  by  analyzing  the  length  of  the  two-component  discrimi- 
nant vectors  were  to  clearly  separate  earthquakes  from  presumed  explosion 
events  from  the  two  regions.  Assuming  normal  statistics,  it  is  expected  that 
90  percent  of  the  presumed  explosions  can  be  detected  with  a probability  of 
false  positive  identification  of  4 x 10"^,  Two  outliers  were  obtained;  one 
from  each  presumed  explosion  region.  This  indicates  the  need  for  more 
discriminant  components  in  the  measured  discriminant  vector,  and  possibly 
more  components  in  the  resultant  discriminant  space  needed  to  identify  all 
explosions  by  region  or  medium  type.  Before  considering  use  of  some  of  the 
more  advanced  capabilities  of  SPEED,  significant  improvements  over  these 
baseline  results  should  be  firmly  established. 
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A brief  summary  of  a series  of  system  studies  describes  the 
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active Seismic  Processing  Module  was  developed  on  a PDP-15  minicomputer 
for  Short-Period  Earthquake/Explosion  Discrimination  (SPEED).  This  pro- 
cessing module  was  imbedded  in  a simulated  special  purpose  seismic  operat- 
ing system.  Interactive  Seismic  Processing  System  (ISPS)  developed  by  Ringdal 
and  Shaub.  All  of  the  capabilities  for  waveform  processing  using  SPEED  were 
described  and  demonstrated  by  processing  a seismic  event  and  displaying  the 
results  of  each  stage  of  processing.  Support  programs  were  developed  to  use 
an  IBM  360/44  computer  to  reduce  all  of  the  station  component  seismograms 
to  a single  representative  event  waveform.  An  option  is  provided  for  inputting 
into  SPEED  the  single  representative  event  waveform  either  corrected  for  or 
uncorrected  for  system  response.  The  examples  analyzed  were  corrected  to 
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The  SPEED  processing  module  as  presently  configured  can  be 
programmed  to  perform  any  one  of  many  discrimination  procedures  by  a pre- 
scribed sequence  of  button  pushing.  To  obtain  baseline  performance  data, 
one  very  simple  and  reproducible  discrimination  processing  procedure  was 
prescribed.  A data  base  of  35  events,  including  20  earthquakes  and  15  pre- 
sumed explosions  in  the  magnitude  range  between  4.4  and  6.  1,  were  run 
through  the  SPEED  processor.  For  each  event,  the  discrimination  measure- 
ments obtained  from  SPEED  were  processed  by  multivariate  discrimination 
analysis  of  frequency  dependent  magnitude  measurements  to  derive  a two- 
component  discriminant  output.  One  component  was  optimized  to  discriminate 
central  Asian  presumed  explosions.  The  other  component  was  optimized  to 
discriminate  western  United  States  presumed  explosions.  The  performance 
results  obtained  by  analyzing  the  length  of  the  two  - component  discriminant 
vectors  were  to  clearly  separate  earthquakes  from  presumed  explosion 
events  from  the  two  regions.  Assuming  normal  statistics,  it  is  expected  that 
90  percent  of  the  presumed  explosions  can  be  detected  with  a probability  of 
false  positive  identification  of  4 x 10  Two  outliers  were  obtained;  one 

from  each  presumed  explosion  region.  This  indicates  the  need  for  more 
discriminant  components  in  the  measured  discriminant  vector,  and  possibly 
more  components  in  the  resultant  discriminant  space  needed  to  identify  all 
explosions  by  region  or  medium  ty’pe.  Before  considering  use  of  some  of  the 
more  advanced  capabilities  of  SPEED,  significant  improvements  over  these 
baseline  results  should  be  firmly  established. 
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SECTION  I 


OVERVIEW 


A.  INTRODUCTION 

Over  the  past  several  years  Texas  Instruments  Incorporated 
has  performed  a series  of  system  studies.  The  first  of  these  defined  the 
surveillance  problem  and  the  functions  required  to  perform  it.  From  this 
it  was  possible  to  describe  facilities,  data  flow,  file  structures,  and  hard- 
ware or  software  requirements  for  seismic  data  processing  and  communica- 
tions processing, 

A report  was  written  to  analyze  the  most  feasible  means  of 
utilizing  past  measurements  of  seismic  events  and  seismic  noise.  Past 
history  of  the  seismic  noise  and  seismic  event  measurements  can  be  used  to 
improve  the  front  end  detection  of  seismic  events  and  to  improve  the  quality 
of  seismic  event  waveform  estimates,  given  the  location  of  a newly  detected 
event.  Another  feasibility  report  was  written  to  evaluate  where  automatic 
data  processing  can  be  most  effectively  used.  Possible  areas  of  application 
are  for  front-end  seismic  event  detectors,  detection  association  processors, 
and  event  waveform  retrieval  processors.  Also  considered  were  the  seismic 
analyst  interpretation  tasks  which  need  to  be  supported  by  interactive  graph- 
ics processing.  This  type  of  information  processing  was  used  for  waveform 
alignment,  later  phase  picking,  recursive  event  focal  parameter  estimates, 
reduction  of  many  event  waveforms  to  one  or  several  to  be  used  for  discrim- 
ination processing,  and  earthquake  and  explosion  discrimination  processing. 

This  report  and  several  preceding  reports  are  part  of  the 
present  phase  of  system  design  studies  to  develop  and  evaluate  the  computer 
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codes  needed  to  perform  surveillance  functions.  To  facilitate  rapid  develop- 
ment, these  are  performed  in  a mode  which  simulates  surveillance  mode 
^ operation.  Interactive  graphics  programs  were  developed  to  perform  long- 

^ period  and  short-period  seismic  event  analysis.  Automatic  programs  were 

developed  and  tested  to  perform  detection  association  processing  and  wave- 
form retrieval. 

B.  SYSTEMS  ANALYSIS  OF  A WORLD-WIDE  SURVEILLANCE 
ENVIRONMENT 

This  report  continues  the  series  of  system  studies  by  Texas 
Instruments  Incorporated  to  develop  feasible  and  economically  viable  options 
of  implementing  a seismic  earthquake  and  explosion  surveillance  system. 

The  first  of  the  reports  by  Sax  et  al.  (1974)  scoped  the  problem  of  defining 
and  parameterizing  a system  capable  of  monitoring  world-wide  occurrences 
of  earthquakes  and  explosions.  One  of  the  parameters  scoped  was  sensor 
deployment,  ranging  from  a network  of  single  sensor  stations  to  a network  of 
array  stations.  Another  parameter  was  the  data  collection  mode.  At  one 
extreme,  a centralized  system  forwards  all  raw  seismic  data  to  a central 
facility.  There,  the  seismic  events  are  detected,  classified  and  documented. 
At  the  opposite  extreme,  to  accomplish  the  same  mission  a decentralized 
system  forwards  only  selected  data  and  information.  In  that  case  data  pro- 
cessing for  event  detection  and  waveform  retrieval  is  partitioned  between 
computer  processing  at  remote  sites  and  computer  processing  at  the  central 
facil  ity . 

The  1974  system  study  report  was  written  to  provide  systems 
analysis  for  implementing  and  operating  a seismic  surveillance  system.  The 
report  specified  definition  of  data  processing  functions,  definition  of  files, 
allocation  of  storage,  control  of  access  to  data  and  information  files  by  auto- 
matic processors  or  analyst-driven  function  processors,  communications, 
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and  communications  control  functions.  Given  this  level  of  function  require- 
ment definition,  it  was  possible  to  specify  preliminary  designs  of  facilities, 
hardware  and  software  modules  needed  to  perform  required  functions,  and 
of  interfaces  needed  to  synchronize  performance  of  the  function  processes. 
The  function  processes  were  designed  to  transform  on-line  files  of  raw 
seismic  data  to  event  associated  data  and  information.  The  function  pro- 
cesses and  file  structures  were  defined  with  sufficient  generality  to  accom- 
modate ground  motion  measurements  either  by  a single  component  sensor, 
multiple  component  sensor,  multiple  sensors  distributed  as  array  elements, 
or  Complexes  of  arra/s.  This  generality  of  software  and  file  structure  en- 
ables the  system  to  adapt  to  changes  in  sensor  deployment.  For  example. 
Single  sensor  stations  could  be  converted  to  array  stations  without  interrrup 
tion  or  serious  adverse  effects  on  the  performance  of  seismic  surveillance. 

If  the  seismic  surveillance  system  is  designed  to  be  imple- 
mented and  operated  as  an  integrated  system,  the  managers  of  the  network 
could  alter  the  size  of  arrays  or  could  add  or  subtract  remote  stations  to 
meet  almost  any  future  anticipated  requirement  for  seismic  event  detection 
capability.  Top  down  planning  of  the  design  for  this  kind  of  flexibility  re- 
quires in  advance  detailed  quantitative  specification  of  all  data  processing 
functions;  of  the  content  and  size  of  storage  files;  of  the  control  needed  for 
sequential  access  to  data  files  by  parallel  operating  function  processes;  of 
the  communications  capacity  requirements;  of  the  system  reliability  and 
recovery  requirements;  and  the  limit  of  time  delay  tolerance  allowable 
before  a final  seismic  event  classification  is  obtained. 

One  important  system  design  decision  is  to  select  either  a 
centralized  or  decentralized  mode  of  seismic  data  collection.  Given  the 
present  state  of  the  art,  either  type  of  data  collection  can  meet  any  reason- 
able anticipated  requirements  for  the  system  monitoring  world-wide  occur- 
rences of  earthquakes  and  explosions.  The  main  impact  of  this  decision  is 
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on  the  size  and  capabilities  of  remote  and  central  facilities  to  hold,  process, 
and  forward  data.  The  decision  should  therefore  be  made  before  firming  the 
design  of  such  facilities.  A decentralized  data  collection  system  utilizes 
station  disk  files  and  remote  station  processors  to  send  reduced  data  and  in- 
formation to  a central  facility.  A centralized  system  utilizes  much  simpler 
direct  sensor  to  satellite  data  collection  links  to  send  all  sensor  data  to  a 
central  facility.  In  both  cases  it  is  anticipated  that  each  remote  station's 
raw  seismic  data  will  be  held  for  at  least  six  to  eight  hours  pending  the  data 
retrieval  needed  to  analyze  each  seistnic  event.  These  event  data  retrievals 
are  driven  by  automatic  and  independently  activated  station  threshold  detec- 
tions. The  detections,  sent  asynchronously  as  bidletins  from  the  remote 
stations  to  a central  facility  arc  linked  by  an  automatic  station  detection 
association  processor  to  specific  seismic  ev'ents  focal  parameters,  e.  g,  , 
origin  time,  location,  and  depth  of  earthquake.  The  detection  association 
processor  utilizes  detection  bulletin  information  to  obtain  the  location  and 
origin  time  of  the  event.  It  must  do  this  to  predict  event  arrival  times  at 
all  remote  stations.  The  predicted  arrival  times  should  be  accurate  to 
within  the  order  of  ten  to  twenty  seconds,  to  assure  confident  retrieval  of 
two-minute  signal  centered  short-period  seismic  event  waveforms  as  well 
as  the  one-half  hour  long-period  waveforms. 

One  advantage  of  a decentralized  over  centralized  data  collec- 
tion is  the  estimated  1 5 to  Z0%  lower  cost  of  implementation  and  10-year  op- 
eration. Another  advantage  is  the  capability  for  .nore  responsive  computer 
assisted  quality  control  of  sensor  operation  at  remote  sites.  The  quality 
control,  e.  g.  , detection  and  repair  of  bad  sensors,  can  be  more  readily  per- 
formed by  staff  personnel  with  equipment  available  at  the  remote  station  pro- 
cessor facilities. 

One  advantage  of  centralized  data  collection  is  simultaneous 
world-wide  access  to  all  of  the  raw  sensor  data.  Because  of  this  access  and 
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the  lack  of  communication  delays,  higher  false  alarm  waveform  retrievals 
can  be  tolerated.  These  would  be  caused  by  detectors  operated  at  lower 
thresholds.  This  could  result  in  some  increased  network  event  detection 
capability.  It  is  also  possible  that  hardware  and  software  economies  of 
scale  can  be  built  into  performing  all  of  the  station  level  detection  and  wave- 
form retrieval  by  means  of  simultaneously  parallel  station  detection  proces- 
sing. Recoverability  from  system  overloading, due  to  abnormal  event  occur- 
rences such  as  earthquake  swarms,  is  faster  and  more  certain  because  of 
simultaneous  access  to  all  of  the  raw  seismic  data  files.  These  can  be  tem- 
porarily dumped  onto  and  retrieved  from  backup  tape  or  disk  units.  If  a 
satellite  data  collection  system  is  built,  it  is  expected  that  considerable  ex- 
cess satellite  data  communications  capacity  would  be  available  for  use  by 
other  agencies  on  a when -available  basis.  For  example,  lower  data  rate 
missions,  such  as  collecting  earthquake  prediction  data,  can  utilize  the  ex- 
cess capacity  of  the  satellite  seismic  data  collection  system  to  forward 
measurements  to  a central  facility.  By  so  sharing  the  cost  of  the  data  collec- 
tion system,  there  could  result  a considerable  communications  cost  savings 
for  participating  agencies.  A disadvantage  of  a centralized  data  collection 
system  is  the  longer  and  more  critical  development  path  with  more  techno- 
logical -isks. 

However  the  decision  is  made  between  centralized  and  decen- 
tralized data  collection,  design  studies  of  algorithms  required  to  perform 
surveillance  system's  functions  can  presently  be  performed.  This  can  be 
done  in  a simulation  mode  prior  to  finalizing  their  design.  By  iterative 
design  and  testing,  certain  critical  functions  such  as  detection  association 
processors,  event  classification  processors,  and  station  detection  processors 
can  be  specified  in  minute  detail.  To  the  extent  that  detailed  specification  of 
function  processors  are  realistic,  the  code  required  to  perform  the  functions 
can  be  developed  free  of  faults.  Finally,  top  down  design  of  an  integrated 
surveillance  system,  i.  e.  , facilities,  hardware,  and  software,  is  feasible 
after  achieving  such  a level  of  problem  definition. 

1-5 


t 


c. 


FEASIBILITY  ANALYSIS  OF  PARAMETER  UPDATE  PROCESSING 


To  optimize  the  operation  of  a network  of  event  detectors, 
thresholds  can  be  controlled  to  obtain  the  highest  possible  detection  capabil- 
ity in  areas  of  surveillance  interest.  For  a network  of  array  stations,  per- 
formance optimization  also  involves  improving  event  estimation  by  correct- 
ing for  consistent  bias  in  beamforming  parameters,  i.  e.  , for  anomalous 
transmission  time  delays,  and  for  anomalous  transmission  amplitude  atten- 
uation factors. 

The  performance  of  the  network  for  acquisition  of  seismic 
event  waveforms  depends  on  the  operating  characteristics  of  the  station  de- 
tectors. These  flag  possible  arrivals  of  events  by  means  of  bulletins.  Given 
these  bulletins,  a detection  association  processor  times  and  locates  the  seis- 
mic events.  At  the  station  level,  sensor  data  is  registered  and  held  on  disk 
pending  the  request  of  event  waveform  data.  A front  end  detection  process 
is  performed  on  the  newly  registered  data.  The  detector  searches  for  indi- 
cations of  new  events  by  detecting  large  amplitudes  on  single  sensors  or  on 
waveform  estimates  derived  from  array  sensors.  The  start  time  of  large 
amplitude  excursions  provides  a first  estimate  of  the  event  arrival  time. 

To  estimate  possible  event  arrivals  on  arrays,  a set  of  targeted  locations 
are  scanned  for  events  originating  from  the  vicinity  of  one  of  the  targeted 
locations.  Thus  array  detections  yield  not  only  event  arrival  times  but  also 
associated  location  ellipses.  In  that  case,  any  point  within  a large  elliptical 
area  centered  about  the  targeted  location  is  a possible  location  of  the  detected 
event.  Since  large  amplitude  excursions  can  occasionally  result  from  ambi- 
ent seismic  noise,  thresholds  are  set  high  enough  to  eliminate  most  large 
amplitude  excursions  due  to  noise.  But  this  is  done  at  the  cost  of  missing 
smaller  events.  Each  level  that  a threshold  can  be  set  can  be  characterized 
by  a rate  of  false  alarms  due  to  noise  as  well  as  by  the  probability  of  missing 
an  event  of  some  specified  magnitude  and  distance.  Each  threshold  level  of 


a station  detector  implies  a point  on  the  curve  of  false  alarm  rate  versus 
probability  of  missing  the  targeted  signal.  This  curve  is  the  operating  char- 
acteristic of  the  front  end  detector.  It  is  important  to  use  thresholds  to  limit 
the  expected  number  of  false  alarm  bulletins.  At  the  network  level,  arrival 
times  at  four  or  more  stations  are  needed  to  time  and  locate  events.  This 
focal  estimate  is  used  to  drive  requests  for  event  waveforms.  This  sampled 
event  data  is  obtained  from  the  raw  seismic  data  held  at  each  station.  The 
false  alarm  bulletins  can  interfere  with  this  process.  By  controlling  the 
average  rate  of  false  alarms  from  all  of  the  stations,  the  network  level  tim- 
ing and  location  of  events  by  detection  association  processing  can  be  optimized. 

Unger  (1974)  analyzed  the  methodology  of  optimizing  threshold 
settings  to  improve  the  detectability  of  events  by  a network  of  sensors.  This 
was  done  by  setting  up  a cost  function  to  assess  the  cost  in  network  perform- 
ance of  detecting  false  alarms  and  of  missing  signals.  The  optimum  threshold 
setting  at  a station  is  that  setting  which  minimizes  the  total  cost  of  both  types 
of  detection  errors.  In  general,  the  false  alarm  rate  of  these  minimum  cost 
threshold  settings  will  not  be  the  same  at  every  station.  For  example,  one 
means  of  implementing  the  minimum  cost  threshold  strategy  is  to  set  up  tar- 
get locations  in  areas  of  surveillance  interest.  In  that  case,  the  false  alarm 
rate  of  the  optimum  threshold  settings  will  depend  on  the  following  parameters; 
the  location  of  targeted  regions,  the  minimum  desired  magnitude  event  to  be 
detected  from  each  region,  and  the  maximum  number  of  false  alarms  that  can 
be  tolerated  to  avoid  missing  events  from  the  region.  Simulation  of  network 
operations  should  allow  these  network  level  threshold  control  parameters  to 
be  set  at  values  which  optimize  the  performance  of  the  detection  association 
processor  in  detecting  seismic  evpnts  from  each  region.  The  detection  bulle- 
tins would  be  coded  to  indicate  each  conditional  region  location  implied  by  the 
threshold.  An  obvious  result  of  the  strategy  of  making  thresholds  conditional 
to  event  location  is  that  thresholds  will  be  set  high  enough  to  prevent  needlessly 
high  false  alarm  rates  due  to  local  events.  This  follows  from  the  extremely 
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high  likelihood  of  detecting  the  targeted  event.  Thresholds  will  also  be  set 
very  high  on  far  out  noisy  stations  due  to  their  very  low  likelihood  of  detec- 
tion. Preliminary  analysis  of  a minimum  cost  threshold  strategy  compared 
to  operating  thresholds  at  a constant  false  alarm  rate  indicated  that  the  mini- 
mum cost  mode  operating  at  the  same  average  false  alarm  rate  can  gain  0.  3 
units  of  magnitude  detection  capability. 

The  operating  characteristics  of  station  level  array  detectors 
depend  not  only  on  threshold  control  but  also  on  the  signal -to-noise  (S/N) 
ratio  of  events  estimated  by  beamforming.  The  S/N  ratio  of  beams  depends 
on  accurate  time  alignment,  the  relative  effect  of  transmission  attenuation  on 
each  sensor,  and  the  relative  noise  level  of  each  sensor.  The  same  noise 
information  used  to  control  thresholds  could  be  used  to  estimate  sensor  noise. 
The  information  needed  to  estimate  the  amplitude  attenuation  and  time  delays 
depends  strongly  on  the  location  and  depth  of  the  event  within  the  targeted  re- 
gion. By  analyzing  past  events  from  each  targeted  region,  systematic  correc- 
tions of  plane  wave  models  for  the  transmission  effects  (c.  g.  , amplitudes, 
velocity,  direction,  curvature,  random  delays,  etc.)  can  be  derived  for  each 
targeted  region.  Such  corrections  have  been  shown  to  be  very  effective  at  the 
Large  Aperture  Seismic  Array  (LASA)  and  the  Norwegian  Seismic  Array  (NOR- 
SAR)  for  estimating  events  of  known  location.  It  is  not  as  yet  known  whether 
such  corrections  can  significantly  improve  the  waveform  estimations  of  knovTi 
events  on  small  arrays  or  significantly  improve  the  performance  of  front  end 
small  array  detectors. 

D.  FEASIBILITY  ANALYSIS  OF  AUTOMATIC  AND  SEISMIC  ANALYST 
INTERACTIVE  DATA  PROCESSING 

A study  was  performed  by  Sax  (1974)  to  evaluate  how  automatic 
and  analyst  interactive  processing  is  most  effectively  employed  for  surveil- 
lance operations.  At  the  present  state  of  the  art,  the  detection  association 
process  severely  limits  the  event  detection  capability  of  a network.  The 
process  is  driven  by  front  end  detection  bulletins  sent  to  a central  facility. 
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There,  false  alarm  bulletins,  large  timing  errors,  and  interfering  transmis- 
sions must  be  eliminated  effectively  from  the  bulletins.  The  remaining  detec- 
tion bulletins  can  then  be  used  to  reliably  detect  and  update  event  focal  param- 
eters. These  updated  estimates  of  the  location,  depth,  and  origin  time  are 
used  to  estimate  expected  arrival  times  of  an  event  at  each  station.  The 
arrival  times  must  be  estimated  to  within  about  one  quarter  minute  to  confi- 
dently retrieve  the  event's  waveform  from  a station.  To  attain  a low  missed 
signal  probability,  front  end  detectors  must  be  operated  at  the  lowest  possi- 
ble thresholds.  To  achieve  this,  the  detection  association  processor  must 
efficiently  eliminate  a large  number  of  false  alarms.  It  is  anticipated  that,  to 
attain  the  most  efficient  elimination  of  false  alarms,  the  detection  association 
processor  will  need  to  be  operated  in  the  automatic  mode  with  a sufficiently 
well  designed  algorithm. 

To  analyze  event  data  after  retrieving  the  waveforms,  it  is 
anticipated  that, due  to  the  interpretive  nature  of  discrimination  processing, 
seismic  analysts  will  be  needed.  These  analysts  will  need  to  turn  over  pro- 
grams efficiently  with  results  displayed  at  each  step.  Interactive  graphics 
programs  are  needed  to  keep  up  with  the  work  required  to  parameterize, docu- 
ment and  classify  events.  At  the  post  detection  association  processing  stage, 
event  waveforms  sent  in  from  different  stations  are  deposited  on  a central 
storage  unit.  The  waveforms  from  each  station  are  then  displayed  for  edit 
by  an  analyst.  The  analyst  selects  waveforms  with  visible  signals  and  time 
aligns  them  by  cursor  manipulation.  This  is  done  interactively  by  the  analyst 
so  that  the  updated  event  location  and  error  ellipse  can  be  observed  after 
timing  each  waveform.  Large  timing  errors  can  be  detected  by  large  changes 
in  the  size  of  the  error  ellipse.  Such  problem  waveforms  can  be  re-timed  or 
deleted  if  large  apparent  timing  errors  are  unresolved.  Another  user  option 
is  to  pick  possible  later  phases,  especially,  nearly  time  coincident  pP  phases. 
As  each  depth-phase  (pP)  is  picked,  the  old  and  new  value  of  the  depth  of 
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focus,  origin  time  estimate,  and  standarfi  deviations  are  displayed.  Converg- 
ence of  the  error  ellipses  can  be  used  for  much  more  accurate  determination 

of  the  focal  parameters,  especially  depth  of  focus. 

In 

After  accurately  determining  the  focal  parameters  of  the  event, 
the  large  S/N  waveforms  can  be  selected  and  corrected  for  the  ground  dis- 
placement system  response,  and  for  absorption.  Each  waveform  can  be  clas- 
sified as  positive  or  negative  first  motion,  with  the  results  stereo  plotted  to 
estimate  the  source  mechanism. 

The  final  step  of  the  post  detection  association  process  is  for 
the  seismic  analyst  to  reduce  the  waveforms  to  one  or  several  to  be  analyzed 
by  the  event  classification  analyst.  The  event  classification  analyst  interac- 
tively performs  extensive  discrimination  analysis  to  classify  the  event  as  an 
earthquake,  as  a possible  explosion,  as  a multiple  explosion,  or  as  a possible 
explosion  hidden  by  the  coda  of  an  earthquake.  For  routine  event  classifica- 
tion analysis,  statistical  discriminants  such  as  M versus  m,  and  variable 

s b 

frequency  magnitude  are  measured  in  order  to  detect  possible  explosions. 
These  statistical  discriminants  arc  supplemented  with  complex  cepstrum 
analysis  to  detect  explosion  echo  delays,  the  reflection  coefficient  of  the 
echo,  the  displacement  waveform  with  the  echo  removed,  and  the  log-log 
spectrum  on  which  the  analyst  picks  the  long-period  amplitude  and  corner 
frequency.  Also,  provision  is  needed  to  analyze  multiple  explosions.  It  is 
expected  that  the  routine  event  classification  analysis  will  be  performed  in 
considerably  less  than  one-half  hour  per  event.  Since  there  is  so  much  of 
this  analysis  which  is  also  highly  iiiterpretive,  it  is  anticipated  that  it  must 
be  done  in  the  analyst-computer  interactive  mode  given  our  present  state  of 
the  art. 

Inevitably,  there  will  be  priority  events  from  areas  of  high 
surveillance  interest  which  will  require  extended  event  classification  pro- 
cessing, This  will  take  considerably  more  interpretation  time  and  more 
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data  analysis  than  routine  event  classification  processing.  For  such  cases  a 
regional  file  will  be  opened  to  document  possible  treaty  violations,  and  to 
analyze  and  evaluate  strategically  located  earthquakes  which  may  present 
possible  opportunities  for  hiding  explosion  sources.  It  is  anticipated  that  a 
need  exists  for  the  capability  of  handling  data  collection  for  such  cases  at  a 
rate  of  approximately  one  per  day.  This  task  would  require  that  a research 
seismologist  interactively  retrieve  the  interesting  event  waveforms  and  other 
associated  nearby  events  for  comparative  analysis.  An  important  support 
task  for  analyzing  interesting  events  is  to  collect  raw  sensor  data  of  consid- 
erable time  duration  from  backup  tapes.  These  extended  waveform  files 
should  be  available  for  deposition  in  the  data  bank  as  interesting  event  files 
within  something  like  twenty-four  to  forty-eight  hours  after  the  event  occur- 
rer  ce.  It  is  reasonable  that  data  supporting  interesting  events  will  be  held 
for  a specified  period  of  time  on  special  event  tapes.  These  would  be  batch 
processed  until  a final  classification  is  obtained.  It  is  anticipated  that  back- 
up tapes  holding  the  raw  sensor  data  will  be  held  at  the  station  levels  for  at 
least  several  months  before  recycling. 

E.  SYSTEM  DESIGN  ANALYSIS  BY  SIMULATION 

To  be  able  to  perform  functions  in  the  surveillance  mode,  it 
is  necessary  to  develop  fault  free  codes  for  all  of  the  automatic  and  interac- 
tive processes  involved  in  detecting  events.  The  capability  to  perform  a 
function,  in  the  surveillance  mode,  can  be  evaluated  by  realistically  simulat- 
ing the  workload  handled  by  each  process.  When  feasible,  this  is  done  by 
analyzing  seismic  event  data  and  seismic  noise  data.  If  such  data  is  not 
available  as  in  the  case  of  hypothetical  fuhire  network  configurations,  it  is 
necessary  to  realistically  simulate  earth  seismicity  and  noise  as  well  as  the 
system  function  to  be  evaluated.  This  is  done  by'  means  of  computer  generated 
Monte  Carlo  realizations  from  representative  statistical  distributions. 
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Statistical  distributions  functions  and,  in  some  cases,  simula- 
tion models  were  used  to  generate  the  effects  of  event  occurrences,  transmis- 
sion, and  the  operation  of  automatic  system  functions  such  as  station  detectors 
and  detection  association  processors.  Shoup  and  Sax  (1974)  constructed  simula- 
tors to  model  the  earthquake  seismicity  and  transmission  of  P waves,  pP,  and 
several  other  dominant  compr es sional  wave  phases.  Also  an  algorithm  was 
used  to  generate  the  statistical  fluctuations  due  to  the  coda  of  earthquakes. 

The  model  used  to  generate  coda  was  based  on  observations  of  the  decay  char- 
acteristics of  earthquakes.  Station  detectors  were  simulated  by  controlling 
the  average  false  alarm  rate  of  each  detector  level,  e.  g.  , of  one  false  alarm 
per  hour.  The  required  threshold  to  attain  the  prescribed  false  alarm  rate 
is  an  output  of  the  station  detector  simulator.  The  event  timing  errors  were 
simulated  by  modeling  a special  type  of  event  complexity  factor:  the  time 

delay  from  first  motion  to  the  maximuna  signal  peak.  Timing  errors  were 
generated  by  automatically  simulating  noise  detection  by  backing  up  from  the 
signal  peak  to  what  was  detected  as  noise  preceding  the  signal.  From  all  of 
the  world-wide  stations  the  detection  bulletins  due  to  false  alarms,  P phases, 
later  phases  and  coda  were  printed  and  copied  on  tape.  From  this  informa- 
tion, the  estimated  performance  and  precision  of  the  simulated  power  detec- 
tors could  be  assessed. 

In  order  to  develop  suitable  codes  for  automatic  detection  as- 
sociation processing,  the  above  tape  containing  world-wide  station  detection 
bulletins  was  used  to  drive  a baseline  model  detection  association  processor. 
The  baseline  procedure  used  was  to  detect  overlapping  location  error  ellipses, 
update  the  location,  and  perform  logical  tests  to  eliminate  false  and  redundant 
information.  The  results  indicated  that,  although  the  bulletin  information 
could  be  used  to  locate  events  and  retrieve  data,  serious  problems  were 
encountered  by  overloading  buffers  and  occasionally  committing  logical  er- 
rors in  the  association  process.  Performance  was  thereby  degraded  by  0.  3 
magnitude  units  of  detection  capability  and  by  the  high  rate  of  about  twenty- 
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five  erroneous  waveform  retrievals  per  day.  Snell  and  Sax  (1976)  developed  a 
more  advanced  model  of  the  same  automatic  detection  association  processor. 
This  was  accomplished  through  repeated  simulations  and  diagnoses  of  results. 
Performance  evaluation  of  detection  association  was  done  by  comparing  de- 
tection association  process  locations  to  the  known  location  of  simulated  events, 
for  all  events  which  exceeded  threshold  at  four  or  more  remote  stations.  The 
new  model  detection  association  processor  demonstrated  almost  no  loss  of 
magnitude  detection  capability.  This  result  was  obtained  at  an  error  rate  of 
about  three  erroneous  waveform  retrievals  per  day  and  at  constant  front  end 
detector  false  alarm  bulletin  rate  of  twelve  per  day  per  remote  station.  The 
original  baseline  detection  association  processor  was  sliown  to  perform  opti- 
mally with  remote  stations  operating  at  the  twelve  per  day  bulletin  false 
alarm  rate.  Since  the  new  model  performs  so  well,  the  false  alarm  rate  of 
front  end  detectors  can  obviously  be  raised  substantially.  This  will  improve 
the  overall  system  detection  capability.  The  optimum  value  of  this  rate  has 
not  as  yet  been  determined  by  simulation. 

After  retrieving  event  waveforms,  the  rate  at  which  the  event 
data  must  be  routinely  interpreted  is  expected  to  range  from  one  to  five  events 
per  hour.  To  keep  up  with  this  flow  of  event  data,  seismic  interpreters  must 
process  between  two  and  four  events  per  hour.  To  keep  the  maximum  event 
processing  delays  reasonably  low,  event  processors  should  be  designed  to 
achieve  rates  close  to  the  upper  limit  but  with  options  for  more  thorough  pro- 
cessing if  the  workload  will  permit  it.  It  is  also  anticipated  that  for  recovery 
purposes,  at  least  one  backup  graphics  processor  will  be  available  to  work  off 
any  occasional  large  backlog  of  event  data.  During  normal  periods,  this  extra 
graphics  processor  would  be  available  to  create  special  data  files  and  to  per- 
form extended  event  classification  analysis  of  the  interesting  events. 

To  normally  maintain  small  backlogs  and  to  otherwise  maintain 
recoverability  from  system  overloading,  it  is  necessary  to  provdde  the  seismic 
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analysts  with  computer  graphics  capability  to  interpret  the  seismic  data.  It 
must  be  possible  for  the  seismic  data  analyst  to  create  required  graphics 
supported  interactive  procedures  by  means  of  a convenient  special  seismic 
command  language.  If  the  graphics  processes  are  properly  designed,  the 
analyst  can  efficiently  select  data  and  processing  tasks,  input  parameters, 
perform  time  series  measurements,  and  recover  conveniently  from  his  own 
interpretation  errors.  The  proper  design  of  the  seismic  graphics  processor 
needed  to  perform  any  analyst-driven  surveillance  function  requires  thorough 
understanding  of  the  surveillance  problem,  the  associated  seismology,  and 
the  efficient  numerical  algorithms  needed  for  the  time  series  analysis.  These 
must  be  integrated  into  a top  down  design  by  the  principal  scientist  familiar 
with  the  above  requirements  working  in  close  conjunction  with  one  or  more 
persons  who  perform  the  system  analysis,  the  computer  science,  and  the 
software  programming.  From  planning  through  to  final  implementation  and 
demonstration,  it  is  usually  necessary  that  these  skills  remain  integrated 
and  coordinated  as  a team  effort. 

In  operating  a surveillance  system  it  is  anticipated  that  at 
least  three  seismic  analysts  arc  needed  in  the  loop  required  to  carry  through 
data  acquisition,  interpretation,  and  deposition  of  a documented  event  report 
into  the  seismic  data  bank.  Initially,  event  waveforms  are  held  on  a central 
storage  element  after  automatic  retrieval  by  the  detection  association  pro- 
cessor, the  station  processors,  and  the  communication  processors.  The 
post-detection  association  seismic  analyst  selects  usable  waveform  data, 
aligns  the  waveforms,  picks  later  phases  such  as  pP,  and  accurately  deter- 
mines focal  information.  As  a final  step,  the  analyst  reduces  the  aligned 
waveforms  to  one  or  several  waveforms  to  represent  source  ground  motion. 

At  this  point  the  routine  event  classification  seismic  analyst  identifies  the 
event  as  an  earthquake,  a possible  explosion,  or  an  interesting  event  requir- 
ing further  research  investigation.  A research  support  seismic  analyst  per- 
forms quality  control  on  the  resultant  event  waveforms  and  information  and 
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deposits  the  classified  events  into  a seismic  data  bank.  This  analyst  supports 
events  of  special  interest  by  retrieving  additional  sensor  data  and  associated 
event  data.  Separate  files  arc  set  up  for  events  of  special  interest,  and  these 
are  possibly  also  deposited  in  the  seismic  data  bank. 

To  establish  some  of  the  interactive  graphics  capability  needed 
to  support  the  surveillance  functions  of  seismic  data  analysts,  Ringdal  and 
Shaub  (1974)  developed  and  demonstrated  an  interactive  long-period  seismic 
processing  system.  They  developed  a high  level  special  purpose  seismic 
command  language  for  fetching  seismic  data  and  performing  any  desired  se- 
quence of  interactive  computer  programs  on  the  data.  Their  system  can  be 
used  to  create  data  processing  procedures  by  performing  serial  or  branched 
sequences  of  commands.  Each  command  allows  the  seismic  analyst  to  ref- 
erence video  displays  of  data  and  information  and  interact  through  a keyboard 
to  guide  the  next  processing  step.  By  selecting  appropriate  processing  steps, 
the  seismic  data  is  manipulated  to  perform  filtering,  magnitude  measurement, 
and  group  velocity  analysis.  These  are  some  of  the  useful  tools  which  can  be 
used  by  the  seismic  analyst  to  perform  routine  event  classification  process- 
ing. In  particular,  filtering  improves  detection  of  the  long-period  event  wave- 
forms. By  measuring  the  long- period  magnitude  of  Rayleigh  waves  and  Love 
waves,  most  explosions  can  be  separated  from  earthquakes.  The  two  random 
populations  are  separated  by  about  0.  7 to  1.0  magnitude  units  with  standard 
deviations  of  from  0.  3 to  0,  45  magnitude  units.  The  detectability  is  approxi- 
mately 1. 0/0.  33*  3.  0.  Given  the  detectability  of  about  3,  at  the  90%  level  of 
explosion  detection,  it  is  expected  that  approximately  5%  of  the  earthquake 

population  will  be  erroneously  classified  as  explosions  based  on  M versus  m 

s b 

criteria  alone.  At  the  50%  level  of  explosion  detection,  the  error  rate  uill  be 
about  one  per  thousand.  Clearly,  better  operating  characteristics  than  these 
are  needed  to  effectively  monitor  a test  ban  treaty.  The  dependence  on  M 

s 

versus  m^  criteria  alone  is  further  complicated  by  repeated  observations  of 
explosion- like  earthquakes  from  certain  regions.  Not  only  is  this  criteria 
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region  dependent  but  at  teleseismic  distances  it  is  only  reliable  for  large  events 
of  the  order  of  5.  0 and  greater  due  to  the  severe  limitations  imposed  Oy  long- 
period  seismic  noise  on  the  detectaljility  of  explosion  surface  waves.  The  a- 
bove  cited  operating  characteristics  are  conditional  on  detecting  both  the  sur- 
face waves  and  the  bodywaves  of  both  classes  of  events. 

Our  present  research  is  motivated  by  the  need  to  effectively 
extend  analyst  interactive  discrimination  capability  to  lower  magnitude  events, 
e.  g.  , on  the  order  of  4.0  magnitude  at  teleseismic  distance,  and  the  need  to 
more  reliably  discriminate  explosion  events  from  earthquakes.  Bache  et  al. 
(1975)  thoroughly  evaluated  a short-period  variable  frequency  magnitude  dis- 
criminant described  by  Savino  and  Archambeau  (1974),  Bache  et  al.  (1974), 
and  Archambeau  et  al.  (1974).  The  results  of  the  evaluation  of  the  short- 
period  variable  frequency  magnitude  discrimination  techniques  indicated  that 
approximately  the  same  operating  characteristics  as  M versus  m^  discrimi- 
nation could  be  extended  to  lower  magnitude  events. 

Any  statistical  method  of  discriminating  earthquakes  and  ex- 
plosions such  as  M versus  m,  measurements  or  short-period  variable  fre- 
^ s b 

quency  magnitude  measurements  will  generate  out  of  the  earthquake  popula- 
tion occasional  and  possibly,  from  some  locations,  consistent  false  indica- 
tions that  an  explosion  has  occurred.  Another  even  more  serious  type  of  er- 
ror is  discriminants  generated  out  of  the  presumed  explosion  population  which 
are  consistently  misidentified  as  ea rtliquakes.  In  order  for  the  surveillance 
system  to  be  credible,  the  rate  of  such  alarms  must  be  adequately  controlled 
j as  part  of  the  discrimination  process. 

\ 

f In  order  to  improve  the  operating  characteristics  of  statistical 

L discriminants,  the  statistical  methods  should  be  backed  up  by  a capability  to 

f separate  the  direct  transmission  and  free  surface  echo.  This  provides  the 

I 

j-  means  of  directly  verifying  that  an  event,  statistically  indicated  as  an 


explosion,  satisfies  the  necessary  conditions  of  an  explosion,  i,  e.  , an  echo 
reflection  coefficient  close  to  minus  one,  a reasonable  time  delay  for  a 
shallow  focus  contained  explosion,  and  a reasonable  direct  transmission 
waveform  shape  and  spectrum. 
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DESIGN  SPECIFICATIONS  OF  THE  SHORT-PERIOD  EARTHQUAKE/ 
EXPLOSION  DISCRIMINATOR  (SPEED)  AS  AN 
INTERACTIVE  GRAPHICS  EVENT  CLASSIFICA TION  MODULE 

ON  THE  PDP-I5 


A.  INTRODUCTION 

The  Short-Period  Ea rthquake/ Explosion  Discriminator  (SPEED) 
processing  module  has  been  integrated  into  the  Interactive  Seismic  Processing 
System  (ISPS)  graphics  capability  on  the  PDP-15  computer.  This  section  pro- 
vides background  information,  data  preparation  procedures,  and  a description 
of  the  SPEED  module. 

In  Subsection  B,  the  ISPS  is  briefly  discussed  in  terms  of  a 
world-wide  surveillance  system.  The  data  preparation  procedures  - the 
reduction  to  one  or  at  most  several  event  waveform  representations  to  be 
used  in  the  discrimination  processing  - are  described  in  Subsection  C.  A 
thorough  description  of  the  operation  of  the  SPEED  module,  with  its  real 
cepstrum  and  Variable  Frequency  Magnitude  (VFM)  and  corner  frequency 
analysis  procedures  is  given  in  Subsection  D.  A general  summary  of  the 
algorithms  used  by  SPEED  are  shown  in  Subsection  E.  Finally,  hardcopy 
examples  of  the  SPEED  graphics  displays  are  shown  in  Subsection  F. 

B,  CONTEXT  OF  IMPLEMENTING  A WORLD-WIDE  SEISMIC 
SURVEILLANCE  SYSTEM 

The  overall  objective  of  a seismic  surveillance  network  is  the 
detection,  location,  and  identification  of  seismic  events  on  a world-wide 
basis.  To  achieve  this  objective,  there  are  six  broad  functional  categories 
that  need  to  be  considered.  In  Table  II- 1,  each  of  the  broad  functional  areas 
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SURVEILLANCE  NETWORK  FUNCTIONS 
{PAGE  1 OF  2) 


Functional  Category 


Information  Generation  and  Processing 


Data  Collection 


Deploy  remote  sensors. 

Data  collection  processing. 

Hold  sensor  data  for  retrieval. 

Edit  event  waveforms  from  the  sensor 
data  which  are  held  for  retrieval. 


Station  Event  Detection 


• Generate  station  detection  bulletins. 


Communications  and 
Data  Lines 


Retrieve  raw  data  from  sensors  to  disk 
or  station  collection  platforms. 
Detection  bulletin  frcim  data  collection 
platform  to  central  storage  element. 
Event  waveforms  from  station  storage 
element  to  central  storage  element. 
Event  parameters  and  representative 
waveform  from  central  storage  element 
to  data  bank. 

Backup  data  tapes. 

Messages  for  command  and  c ontrol  of 
central  and  station  operations. 
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TABLE  II-l 

SURVEILLANCE  NETWORK  EUNCTIONS 
(PAGE  2 OE  2) 


E'unctional  Category 


Retrieval  of  Event 
W aveforms 


Storage  and  Retrieval 
of  Event  Parameters 
and  Data 


System  Control  and 
Upgrade 


Information  Generation  and  Processing 


Control  retrieval  of  event  data  from 
station  storage  elements  by  generating 
event  locations  and  origin  times  from 
detection  bulletin  information. 
Generate  accurate  focal  parameters 
from  event  data  sent  in  from  stations. 
Classify  events  by  discrimination 
analysis  of  representative  event  wave- 
forms. 


File  event  bulletins  and  representative 
waveforms. 

File  hourly  noise  parameter  measure- 
ments. 

Update  regional  reference  event  file. 
Update  and  consolidate  index  and  files 
for  long-term  storage  of  events. 


Stations  and  network  performance 
evaluation. 

Prediction  and  control  of  system  over- 
loading. 

' Fault  detection,  maintenance  and  re- 
pair control. 

< Start  up  and  updating  of  new  model 
function  processors. 

I Updating  reference  event  files. 

Research  support  by  data  base  genera- 
tion. 


are  defined.  They  are  broken  down  into  specific  function  processes  which 
generate  the  files  of  information  needed  to  control  the  operation  of  a seismic 
surv'eillance  network. 

The  hardware  configuration  of  a seismic  surveillance  network 
is  shown  in  Figure  II- 1.  Arrays  or  regional  deployments  of  short-period  and 
long-period  sensors  are  linked  by  communications  to  a data  collection  plat- 
form. At  the  data  collection  platform,  the  data  are  formatted,  multiplexed, 
and  held  temporarily  on  a station  disk  and  permanently  on  a backup  tape. 
Station  detection  functions  are  immediately  performed  upon  this  sensor  data 
held  i.n  the  station  disk.  Independently  of  other  stations  in  the  network,  each 
station  detection  processor  detects  the  possible  arrival  of  seismic  events. 
These  detections  are  described  on  detection  bulletins  which  are  immediately 
sent  by  a communication  link  or  data  line  from  the  data  collection  platform 
to  a central  facility.  At  the  central  facility  the  detection  bulletin  is  decoded 
by  the  COM  processor  and  queued  on  the  central  disk  system  in  a detection 
bulletin  file.  On  a first  in/first  out  basis,  the  detection  bulletins  fill  the  work 
space  of  an  automatic  detection  association  processor.  This  processor  uses 
inform.ation  obtained  from  the  incoming  detection  bulletins  to  generate  output 
which  includes  locations  and  origin  times  of  events.  The  event  focal  param- 
eter estimates  are  used  at  this  point  to  automatically  request  time  windows  of 
data  containing  the  event  waveforms.  These  requests  to  retrieve  a seismic 
event  waveform  data  are  sent  by  communication  links  or  data  lines  to  the 
station  data  collection  platforms.  At  each  station  data  collection  platform, 
the  desired  event  waveforms  arc  automatically  retrieved  from  a disk  storage 
element  which  is  holding  the  seismic  data  for  up  to  six  hours.  The  retrieved 
event  waveform  data  are  sent  back  as  data  packets  during  periods  of  time  that 
the  communication  link  or  data  line  is  available  (e.  g.  , between  transmissions 
of  detection  bulletins).  The  waveform  data  are  received  at  the  central  facility 
by  a central  communications  processor  which  files  the  waveforms  on  the  cen- 
tral storage  element  under  supervision  of  a central  disk  system. 


Seismic  Sensors  DAP  Detection  Association  Processor 

Station  Detection  Processor  ECPl  Event  Classification  Processor  (Discrimination) 


SEISMIC  NETWORK  HARDWARE  CONEIGURATION 


r 


I 


I 
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The  event  data  stored  on  the  central  disk  system  consists  of 
many  channels.  To  be  practical,  discrimination  processing  must  be  fast 
enough  to  keep  up  with  the  data  flow  generated  by  the  surveillance  network. 

The  reduction  of  the  multi-sensor  representations  of  an  event  to  one  or  sev- 
eral representations  for  discrimination  processing  is  described  in  the  next 
section. 

C.  DATA  PREPARATION  FOR  DISCRIMINATION  PROCESSING 

The  event  waveform  file  stored  on  the  central  disk  system  con- 
sists of  one  or  more  sensors  or  signal  estimates  from  each  station  receiving 
the  event.  These  must  be  reduced  to  one  or  at  most  several  event  wav'eform 
representations  to  be  used  for  discrimination  processing.  In  order  to  rapidly 
classify  the  event  as  a possible  explosion  or  earthquake,  the  discrimination 
processing  must  be  rapid  enough  to  keep  up  with  data  flow  on  the  order  of  50 
events  per  day.  .Mlowing  a utilization  of  50%  of  computer  capacity,  the  pro- 
cessing time  per  event  for  all  routine  ev^ent  classification  processing  should 
be  on  the  order  of  15  minutes.  Specific  algorithms  should  therefore  consume 
no  more  than  several  minutes  of  the  event  processing  time.  In  order  to  reli- 
ably classify  the  events,  there  should  be  considerable  diversity  of  discrimi- 
nation processing  readily  accessible  to  the  analyst.  In  particular,  short 
echo  time  determinations  by  cepstrum  analysis  should  be  mandatory  in  order 
to  reduce  false  identification  of  earthquakes  to  an  acceptably  low  false  alarm 
rate  for  events  of  magnitude  greater  than  the  minimum  targetted  for  surveil- 
lance, Also,  some  form  of  the  variable  frequency  magnitude  method  should 
be  used,  provided  that  discrimination  capability  comparable  with  M^  versus 
m^,  can  be  obtained  with  better  detectability  of  lower  magnitude  events. 

At  the  station  platform  level,  multi-sensors  such  as  arrays 
arc  reduced  to  a single  station  waveform  representation.  Since  very  accurate 
arrival  time  measurements  are  not  required  at  this  level,  the  reduction 
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process  can  be  carried  out  automatically.  The  need  for  reduction  also  occurs 
at  the  post  DAP  level,  where  station  waveforms  are  accurately  timed  by  an 
analyst  and  precise  focal  parameter  estimates  are  made.  A post  DAP  analyst 
edits  stations  with  seismic  information,  signal  centers  the  event,  and  invokes 
processing  to  reduce  the  station  waveforms  to  a single  event  waveform  repre- 
sentation to  be  queued  for  discrimination  processing  by  the  event  classifica- 
tion analyst. 

A requirement  for  designing  data  reduction  algorithms  for  the 
abov'e  applications  is  to  produce  event  estimates  with  an  undistorted  spectrum 
of  the  event.  Beamforming  is  in  many  cases  inadequate,  as  deviations  from 
the  assumption  of  plane  waves  in  an  infinite  homogeneous  medium  (the  model 
for  beamforming)  produces  severe  attenuation  of  high  frequencies  which  de- 
grades the  data  for  short-period  discrimination  processing.  A method  of 
event  estimation  based  on  multi-channel  complex  cepstrum  estimation  is 
shown  in  Figure  II-2.  The  critical  technical  points  in  the  development  are  to 
prevent  aliasing  of  the  phase  between  channels,  and  to  develop  the  optimum 
method  of  automatically  weighting  channels  in  computing  the  average  log  am- 
plitude and  phase  of  the  representative  waveform.  In  Figure  II-2,  it  is  also 
shown  that  the  waveforms  are  corrected  for  spectral  bias  due  to  additive 
noise,  and  are  corrected  for  system  response  to  produce  an  estimate  of 
broadband  ground  motion  such  as  acceleration,  velocity,  or  displacement. 

As  a final  step,  single  channel  complex  cepstral  techniques  such  as  short- 
pass  filtering  are  used  to  reduce  the  extensive  coda  of  large  complex  events. 

D.  APPLICATION  OF  THE  INTERACTIVE  SEISMIC  PROCESSING 
SYSTEM  TO  DISCRIMINATION  PROCESSING 

The  Interactive  Seismic  Processing  System  (ISPS)  provides  an 
interactive  graphics  capability  on  the  PDP-15  computer  located  at  the  Seismic 
Data  Analysis  Center  for  the  purpose  of  detecting  and  analyzing  seismic  wave- 
forms. It  currently  consists  of  a system  supervisor  and  five  stand-alone 
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DATA  BASE  INPUT  PREPARATION 


Accurately  Time 
Signals,  Pick 
Phases,  and 
Locate  Event 


J EVENTS  ON  TAPE  WITH  1 WAVEEORM  PER  EVENT 


EIGURE  II-Z 

REDUCTION  OE  WAVEEORMS  TO  PARAMETERS 
AND  A SINGLE  WAVEEORM  REPRESENTATION  OE  THE  EVENT 


processing  modules  as  shown  in  Figure  II- 3.  In  particular,  they  are  defined 


as  follows: 

SUPVSR 

The  system  supervisor  to  control  processing  module 
execution  within  the  CRT  console  environment. 

DPSCAN 

- A processing  module  to  catalog  and  display  waveforms 
from  the  large  capacity  disk. 

SELEV 

- A processing  module  to  select  any  event,  station,  com- 
ponent, or  time  window  for  analysis. 

FILTER 

- A processing  module  to  perform  general  signal  analysis 
functions. 

GRVEL 

A processing  module  to  generate  long-period  dispersion 
curves  and  fundamental  mode  source  spectrum. 

SPEED 

- A processing  module  to  perform  short-period  earthquake/ 
explosion  discrimination  utilizing  cepstrum,  variable 
frequency  magnitude,  and  corner  frequency  techniques. 

The  continuing  methodology  associated  with  ISPS  is  the  simula- 
tion of  analysis  functions  to  be  performed  by  processors  comprising  a world- 
wide seismic  surveillance  system.  A description  of  how  this  simulation  is 
achieved  is  given  as  follows; 

• Research  and  development  of  analysis  software  for  detection, 
association,  and  classification  of  seismic  events.  The  compo- 
nents of  this  software  are  computational  subroutines  coded  in 
FORTRAN  IV,  having  minimal  or  non  existent  I/O  require- 
ments. This  insures  that  these  routines  are  machine  independ- 
ent and  are  easily  maintained. 

• Imbedding  this  analysis  software  in  an  environment  simulating 
that  of  a surveillance  system  when  each  processor  is  driven 
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RESEARCH  ON  ROUTINE  EVENT  CLASSIFICATION  PROCESSING 
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FIGURE  II-3 

INTERACTIVE  SEISMIC  PROCESSING  SYSTEM  (ISPS) 
INFORMATION  FLOW 
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by  a special  purpose  operating  system  which  supports  the 
following: 

A console  device  to  serve  as  the  physical  interface  be- 
tween computer  system  and  analyst,  capable  of  commu- 
nicating and  displaying  symbol  strings  and  graphic 
images 

A command  language  to  control  execution  of  analysis 
functions  assigned  to  the  processor 

An  interpretive  high  level  seismic  programming  language 
permitting  the  analyst  to  design  and  execute  special 
functions  on-line.  The  language  will  have  the  capability 
to  access,  perforin  computation  upon,  and  display  data 
for  evaluation  and  analysis  purposes 

Automatic  processing  modes  providing  for  execution  of 
analysis  functions  with  minimal  intervention  required  by 
the  analyst  once  processing  procedures  have  been  stand- 
ardized. 

• Insuring  that  the  ISPS  environment  accurately  simulates  that 

of  an  actual  processor  by  incorporating  a toji-down  design 
stressing  the  following; 

Operational  flexibility 

System  reliability  and  recoverability 

System  maintainability  and  extendability . 

This  effort  was  given  high  priority  for  two  main  reasons.  First,  interactive 
processing  systems  must  be  capable  of  accepting  erroneous  input  without  re- 
quiring a restart  of  the  entire  analysis  procedure  in  question.  Secondly,  the 
system  should  be  extendible  so  as  to  maintain  its  applicability  as  advance- 
ments in  seismic  processing  arc  realized. 


A new  analysis  function  designated  SPEED  has  been  developed 
to  simulate  discrimination  processing  by  utilizing  CEPSTRUM  and  Variable 
Frequency  Magnitude  (VFM)  techniques  and  corner  frequency.  SPEED  will 
be  demonstrated  as  a processing  module  which  has  been  incorporated  into 
ISPS.  Figure  II-4  illustrates  the  utilization  of  a procedure  named  ISPED  to 
perform  discrimination  processing  within  the  ISPS  environment.  ISPED  de- 
fines a processing  module  execution  sequence  which  includes  the  modules 
DPSCAN,  SELEV,  and  SPEED. 

The  function  SPEED  consists  of  three  discriminant  calculations. 
In  F'igure  II-5,  the  first  one  of  the  discriminant  functions  performs  cepstrum 
analysis;  the  second,  variable  frequency  magnitude  (VP'M)  analysis,  and  the 
last,  corner  frequency  analysis  (CF).  The  three  functions  of  SPEED  " cep- 
strum, VFM,  and  CF  - are  described  as  follows: 

• The  cepstrum  analysis  is  subdivided  into  three  interactive  sub- 
routines. SUBl  adjusts  start  time  or  exponentially  tapers  the 
data  as  desired  by  the  analyst.  SUB2  separates  the  signal  and 
echo  by  trial-and-er ror  deconvolution.  SUB3  generates  a re- 
sidual seismogram  and  sets  up  CEPSTRUM  to  pick  later  phases 
(e.g.  , multiple  explosions),  if  desired  by  the  analyst.  The 
statistics  generated  for  discrimination  are  the  time  delay  of 
the  echo,  reflection  coefficient,  correlation  coefficient  be- 
tween the  echo  and  the  first  arrival,  and  the  residual  noise 
when  the  estimate  of  the  echo  and  first  arrival  are  removed 
from  the  original  data, 

• The  VFM  analysis  in  SUB4  removes  the  exponential  taper  from 
the  first  arrival  signal  estimate  derived  from  CEPSTRUM.  In 
SUB5  it  computes  the  high  and  low  frequency  dependent  magni- 
tudes from  the  corresponding  semi -log-signal  spectrum  dis- 
play corrected  for  absorption. 
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The  procedure  ISPED  may  be  created  and  executed 
at  the  Interactive  Console  by  the  following  commands; 


. . . . 

. CREATE  ISP  ED 

1. 0 

. PERFORM  DPSCAN 

2.  0 

. PERFORM INTERUP 

3.  0 

. PERFORM  SELEV 

4,  0 

. PERFORM  INTERUP 

5,  0 

. PERFORM  SPEED 

6,  0 

. PERFORM  INTERUP 

7.  0 

. #S 

. PERFORM  ISPED. 

FIGURE  II-4 

IMBEDDING  SPEED  WITHIN  ISPS  ENVIRONMENT 
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SHORT- PERIOD  EARTHQUAKE/ EXPLOSION  DISCRLMINATOR 


Processing  Moduli. 


SPEI':D 


Sub-Motiuie 

Sulj  - Module 

CEPSTRUM 

VFM 

SUL  1 


SUL  Z ! 1 SUIi  3 


SUL  4 I SUL  : 


SUL  6 


SUL  1 - Conditions  Dita  For  Cep.'^lruni  Analysis 

SUL  Z - Separ;-tcs  Signal  And  Echo 

SUL  3 - Computes  Resitiual  Seismogi'ain 

SUB  4 - Conditions  Data  For  Sp'.clral  Analysis 

SUB  ''  - Computes  Variable  Freque.nry  ^lagnitudes 

SUB  6 - Computes  Corner  Frequency  .And  Amplitude  I'T'orn  Log-Log 

Spectr.irn  Plot 


FIGURE  II-S 

SHORT-PERIOD  EARTHQUAKE/EXPLOSION  DISCRIMINATOR 


• The  CF  analysis  in  SUB6  can  measure  up  to  four  corner  fre- 

quencies, amplitudes,  and  roll-offs  from  the  log  amplitude 
versus  log  frequency  spectrum  of  the  signal  transmission 
corrected  for  absorption. 

The  software  configuration  of  the  interactive  CEPSTRUM  pro- 
cessor is  shown  in  Figure  II-6.  A signal  is  input  to  SUBl,  and  the  data  and 
real  cepstrum  are  displayed.  The  complex  log  spectrum  is  calculated  and 
held  in  workspace  along  with  the  data.  If  the  analyst  is  not  satisfied  with  the 
start  time  of  the  data  selected  with  the  SELEV  module,  he  may  input  other 
trial  time  shifts  until  he  is  satisfied.  The  analyst  then  presses  the  taper  but- 
ton until  he  is  satisfied  with  the  appearance  of  the  data  and  the  real  cepstrum. 

If  he  tapers  too  much,  he  recovers  by  pressing  the  replace  taper  button. 

After  data  conditioning  is  finished,  the  analyst  presses  the  hardcopy-exit 
button  and  a Calcomp  plot  tape  of  the  results  is  prepared.  Should  the  analyst 
desire,  the  data  may  be  plotted  while  he  is  performing  other  analysis.  Fol- 
lowing preparation  of  the  plot  tape,  control  is  passed  to  SUE2.  In  SUB2,  the 
cepstrum  and  conditioned  data  display  are  refreshed.  The  analyst  presses 
the  deconvolve  button.  He  then  inputs  trial  values  of  the  reflection  coefficient 
and  time  delay,  examining  the  effect  of  deconvolution  on  the  data  and  cepstrum. 
When  satisfied  with  the  resultant  deconvolution,  he  presses  the  load  button 
which  computes  and  saves  the  echo  and  allow'S  him  to  short  pass  the  signal. 
Smoothing  of  the  log  spectrum  of  the  signal  estimate  with  a Bartlett  window 
provides  a short  pass  filter.  This  is  used  to  separate  the  earthquake  pulse 
from  the  echo.  If  the  analyst  presses  the  smooth  button  one  too  many  times, 
he  can  recover  the  preceding  state  by  pressing  the  replace  smooth  button. 

By  pressing  the  load  button  again,  the  smoothed  signal  is  saved  and  the  echo 
is  retrieved  to  be  similarly  short  passed.  The  load  button  is  pressed  a third 
time  to  forrn  the  reconstructed  signal  with  its  echo.  Whenever  the  load  but- 
ton is  pressed,  a display  of  the  conditional  data,  the  signal,  the  echo,  and  the 
reconstructed  signal  arc  shown  from  top  to  bottom  as  they  become  available 
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Separates  signal  from  coda  and  noise 

Separates  the  direct  path  and  free  surface  reflection 

Estimates  first  pulse  direct  path  and  echo 

Sets  up  SPEED  to  compute  multiple  events 

Computes  discriminant  statistics 


by  the  above  analyst  processing.  After  finishing  SUEZ  processing,  the  analyst 
presses  the  hardcopy  - exit  button  to  plot  the  results,  and  passes  control  to 
SUB3,  In  SUB3,  the  reconstructed  seismogram  is  subtracted  from  the  condi- 
tioned data  to  produce  a residual  seismogram  with  the  first  phase  signal  and 
echo  removed.  A later  phase  is  picked  by  means  of  trial  shifts  of  the  wave- 
form similar  to  that  done  in  SUBl.  By  pressing  the  SAVE  PARAMETER  but- 
ton, the  residual  signal,  which  was  shifted  to  the  front  of  the  seismogram,  is 
replaced  as  input  to  CEPSTRUM.  In  this  way,  up  to  four  phases  (c.  g.  , multi- 
ple explosions)  can  be  successively  analyzed  by  the  analyst.  In  the  preceding 
operations,  discriminants  such  as  echo  time,  reflection  coefficient,  correla- 
tion coefficient,  and  residual  signal -to -noise  ratio  are  calculated  and  included 
in  the  hardcopy  as  annotations. 

The  software  configuration  for  variable  frequency  magnitude 
(VFM)  calculation  is  shown  in  Figure  II-7.  Initially  in  SUB4,  the  smoothed 
and  tapered  signal  estimate  from  CEPSTRUM  is  displayed  along  with  the  real 
cepstrum.  The  data  are  conditioned  by  inverting  the  exponential  tapers  put  in 
by  the  analyst  tc  condition  the  data  for  the  cepstrum  analysis.  Roundoff 
».  rrors  which  may  occur  at  the  end  of  the  record  are  zeroed  out  by  pressing 
the  zero  button.  The  signal  may  be  short  pass  filtered  by  pressing  the  smooth 
button.  Control  is  then  passed  to  SUB5  and  the  real  log  spectrum  of  the  sig- 
nal estimate  is  displayed.  This  spectrum  may  be  smoothed  incrementally  by 
a Bartlett  window  and  may  be  corrected  for  absorption.  Magnitude  measure- 
ments are  initiated  by  pressing  the  measure  magnitude  button.  The  average 
log  amplitude  is  computed  in  a central  window  of  the  log  spectrum.  The 
change  in  log  amplitude  between  a marked  low  frequency  and  the  average  and 
a marked  high  frt'quency  and  the  average  are  computed  as  VFM  discriminants. 

If  the  signal  is  obviously  noisy  at  the  specified  low  or  high  frequency,  a cursor 
is  used  to  pick  points  on  the  spectrum  near  the  desired  frequency  outside  the 
frequency  range  of  the  noise.  By  linear  extrapolation,  the  magnitude  difference 


• Operates  on  signal  estimate  output  of  cepstrum  program 

• Inverts  the  taper  which  was  applied  to  the  signal 

• Add  zero's  to  end  of  signal  model,  if  needed 

• Picks  two  low  and  high  frequencies  to  read  magnitude 

• Linearly  extrapolates  to  obtain  low  and  high  frequency 
magnitude  at  desired  frequencies 


SUP  6 rOMPL'Tf  COP  VI  « PHIOUENCY 


FIGURE  II-7 

INTERACTIVE  PROCESSOR  SPEED  MODULES 
VARIABLE  EREQUENCY  MAGNITUDE  AND  CORNER  FREQUENCY 
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at  the  desired  high  or  low  frequency  is  computed.  Our  experience  indicates 
that  this  option  is  usually  not  needed. 

The  software  configuration  of  the  corner  frequency  (CF)  calcu- 
lation is  also  shown  in  Figure  II-7.  The  real  log  spectrum  versus  frequency 
computed  in  SUBS  is  transformed  into  a log  spectrum  versus  log  frequency 
display.  Corner  frequencies,  magnitudes,  and  roll-offs  are  measured  by 
pressing  the  corner  frequency  button.  A cursor  is  used  to  pick  two  points  of 
the  spectrum  on  each  side  of  the  desired  corner  frequency.  The  corner  fre- 
quency is  then  automatically  picked  as  the  intersection  of  the  two  straight 
lines  through  the  left  and  right  pair  of  points.  This  process  may  be  repeated 
up  to  four  times.  The  resultant  corner  frequency,  amplitude  at  the  corner 
frequency,  and  roll-offs  on  the  left  and  right  side  of  the  corner  are  displayed 
for  each  corner  picked.  If  the  analyst  is  not  satisfied  with  a picked  corner 
frequency,  it  may  be  erased  by  pressing  the  replace  corner  frequency  button. 
After  the  analyst  finishes  picking  corner  frequencies,  the  hardcopy -exit  but- 
ton can  be  used  to  produce  a Calcomp  plot  tape. 

An  exit  button  is  present  in  all  of  the  subroutine  modules  de- 
scribed above.  This  is  a recovery  procedure  which  re-starts  the  subroutine 
module  in  case  of  irrecoverable  errors  or  failures.  It  is  unlikely  tiiat  this  op- 
tion will  be  needed  because  of  other  internal  provisions  for  error  recovery. 

E.  ALGORITHMS 

1.  Data  Preparation 

a.  Align  and  Display  Data 

A partial  list  of  symbols  used  in  the  following  derivation  of 
algorithms  is  given  here. 
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represents  uniform  time  between  sampled  data  points, 
number  of  time  sampled  points. 

sequence  number  of  real  time  domain  data  points, 
sequence  number  of  complex  frequency  domain  data  points, 
uniform  frequency  interval  between  spectrum  samples, 

1 /NT. 

fast  Fourier  Transform  operator;  FFT~^(  ),  inverse. 

data;  d(nT),  sampled  data. 

spectrum  computed  from  sampled  data, 

real  part  of  Z(kAf). 

imaginary  part  of  Z(kAf). 

amplitude  spectrum 

phase  spectrum 

the  finite  difference  of  phase  as  an  approximation  to  the 
differential  of  the  phase. 

difference  of  the  real  part  of  the  spectrum, 
difference  of  the  imaginary  part  of  the  spectrum, 
exponential  taper  used  to  derive  channel  Wiener  weights, 
et  al. 

sampled  noise  data  on  channel, 

sampled  signal  data  on  channel, 

parameter  to  weight  ambient  noise  in  deriving  channel 
Wiener  weights. 

parameter  to  weight  coda  noise  in  deriving  channel 
Wiener  weights. 

signal  power;  NF  , noise  power;  C.  , coda  power  used  in 
channel  Wiener  weight  estimates. 

the  spectrum  of  a linear  coda  operator  for  the  channel, 
the  spectrum  of  signal  common  to  all  channels. 
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Nj(f) 

Log  Zg(kAf) 
Log  Ag(kAf) 
(f>g(kAf) 

Log  Ag(kAf) 

c/)|(kAf) 

S(kAf) 

d-(nT) 

R(Z) 

Rj(Z) 

R2(Z) 

a. 

J 

"j 

d'(nT),  d(nT) 
Log  A”(kAf) 
(f)^(kAf) 

Re(  ) 

Log  I \ 
()-►() 

DISPLAY(  ) 
T/Q 


xL, 

the  spectra  of  independent  noise  added  to  the  j channel 

log- spectrum  of  signal  corrected  for  noise  bias. 

real  part  of  log- spectrum  of  signal. 

imaginary  part  of  log- spectrum  of  signal. 

corrected  for  absorption. 

corrected  for  absorption. 

exponentiated  spectrum  corrected  for  absorption, 
estimate  of  signal  pulse  corrected  for  absorption, 
system  response 

part  of  system  response  due  to  poles, 
part  of  system  response  due  to  zeros, 
roots  of  poles  of  system  response. 

scale  factors  of  system  response  associated  with  poles, 
corrected  for  poles  of  system  response, 
binomial  smoothing  of  Log  Ag(KAf). 
binomial  smoothing  of  c/)g(kAf). 
real  part  of  complex  elements  of  an  array, 
logarithm  of  complex  eleinents  of  an  array, 
the  elements  generated  by  operators  on  the  left  side  re- 
place the  elements  of  the  array  on  the  right  side, 
display  of  arrays  designated  in  parenthesis, 
travel  time  divided  by  absorption  Q when  used  in  Strick' 
absorption  formula. 


Given  the  ray  parameter  dt/dr  {seconds/kilometer)  and 
direction  0 of  the  incident  wave. 


dt/dx  = dt/dr  cos  0 


dt/dy 


dt/dr  sin  0 


(II- 1) 


and  the  travel  time  delay  from  the  center  of  the  array  (x  , y ) to  a sensor 

•o  o 

position  (x.,  y.)  is;  i- 

At.  = (dt/dx)  (x.  - X ) + (dt/dy)  (y.  - y ) . (H-2) 

1 1 o 1 o 

A matrix  plot  of  all  of  the  plane  wave  aligned  sensors,  eitlier  filtered  or  unfil- 
tered and  with  the  signal  centered,  is  inspected  by  the  analyst  to  remove  in- 
operative channels,  noisy  channels,  and  clipped  channels, 

b.  Selected  Channels  Are  Complex  Cepstrum  Beamformed 

Acceptable  data  channels  are  timed  to  allow  about  one-half 
second  of  noise  preceding  the  earliest  apparent  arrival  of  signal  energy. 

From  that  start  time,  the  following  512  points  (20  points  / second ) of  data  are 
read  into  a signal  buffer;  the  preceding  512  points  are  read  into  a noise  buffer. 
/\n  exponential  taper,  a'^  - where  n is  the  number  of  points  delay  from  the 

start  point  - is  applied  to  the  data  in  the  signal  and  noise  buffer.  The  param- 
eter, a , weights  the  cepstrum  analysis  more  heavily  toward  the  start  time  of 
the  signal.  A recommended  value  is  0.98  which  is  down  to  0.  1 3 after  100 
points.  If  desired,  snialler  values  can  be  applied  to  the  noise  buffer  for  more 
conservative  estimates  of  the  interfering  noise  power.  We  used  0.  98  for  both 
buffers.  The  signal  buffer  and  noise  buffer  were  filled  out  with  1536  zeros  so 
that  each  buffer  contained  2048  points. 

An  N = 2048  point  fast  Fourier  Transform  (FFT)  was  applied 
to  the  signal  and  noise  buffer  producing  a signal  and  noise  spectrum,  each  of 
1023  complex  points  plus  2 real  points  at  zero  frequency  and  the  Nyquist  fre- 
quency (10  Hz).  The  log- spectrum  is  computed  of  the  spectral  data  in  the 
signal  and  noise  buffer.  Given  the  data,  d(t)  , as  uniformly  sampled  signal  or 
noise  at  period  T,  the  spectrum  sampled  at  frequency  interval  Af  = 1/NT  fol- 
lows : 

Z(kAf)  = FFT(d(nT))  n = l,..,,N. 


11-22 


X(kAf)  + i Y(kAf)  k=0,  ...,N/2 


Log  (Z(kAf))  = Log  A(kAf)  + i</)(kAf) 


(n-3) 


/ 2 2\1 /2 
A(kAf)  = (X(kAf)  + Y(kAf)  ) 


</)(kAf)  = tan"  (Y(kAf  )/X(kAf)) 


FFT{  ) = functional  Fourier  transform  operation. 

The  phase  </)(kAf)  is  complicated  by  rapid  changes  in  phase 
from  one  frequency  point  to  another  due  to  unknown  start  time  delays  and  in- 
fluence of  coda.  This  can  cause  errors  in  detecting  the  position  and  size  of 
2n  jumps  in  the  phase  angle,  which  must  be  corrected  to  produce  a continu- 
ous phase  curve.  To  avoid  this  problemi,  the  phase  information  is  passed 
through  as,  i.  e.  , A<MkAf),  a first  difference  operator  running  from  k=l,  to 
k=N/2,  where  kAf  are  the  corresponding  frequencies  running  from  zero  to 
Nyquist  (N/2)  • (1/NT)  = 1/2T. 


A</>(kAf)  = 


Y(kAf)  AX(kAf)  - X(kAf)  AY(kAf) 
fx(kAf)^  + Y(kAf)^) 


(II-4) 


where 


AX(kAf)  = X(kAf)  - X((k-1  )Af) 


AY(kAf)  = Y(kAf)  - Y((k-1  )Af). 


The  data  on  each  channel  d(nT)  are  converted  to  a log- spectrum, 
where  the  first  part  is  A(kAf)  and  the  second  part  is  Af/>(kAf);  and  both  parts 
are  even  symmetric  at  zero  and  the  Nyquist  frequency  (k  = 0,  N/2). 

Each  channel  ^A(kAf),  A‘/»(kAf)^  is  weighted  and  summed  by 
a Wiener  weight  for  the  channel: 

W.  = S^/(S^  + uN^  + vC^)  (II- 5) 

J J J J J 
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where 


n=  0 ^ 

N 

N-  = E [j-^N'"^’  • ""? 

■’  n=  0 

■'  n=  0 

u = parameter  between  zero  and  one  to  weight  ambient 
noise 

V = parameter  to  weight  coda  measure 
w = exponential  taper,  zero  to  less  than  one. 

For  the  channel  Wiener  weights  to  be  realistic  signal  and  noise,  it  is  recom- 
mended that  w be  set  to  0.  94  to  0.  98.  A value  of  w = 0.  95  was  used  for 
this  study,  which  effectively  gives  a signal  window  length  of  about  two  sec- 
onds. Note  that  setting  u = v = 0 pi'ovides  a unit>^  weighted  beamform  pro- 
cess. Setting  u = 1,  v = 0 provides  a full  Wiener  weighting  using  ambient 

2 2 2 

noise.  It  is  recommended  that,  since  S , N.  and  C.  are  subject  to  statis- 

J 1 J 

tical  errors  as  well  as  bias  due  to  timing  problems,  less  than  full  Wiener 
weights  be  used.  In  this  study,  u = v = 0 in  order  to  baseline  our  results 
using  unity  weighted  complex  cepstrum  beamforming.  It  is  also  noted  that 
the  Ao(kAf)  of  the  noise  buffer  contains  no  information  of  value.  It,  there- 
fore, was  not  computed  and  passed  through  the  next  stage  of  data  processing. 

c.  The  Complex  Cepstrum  Beamform  is  Smoothed  and  Cor- 
rected for  System  Response,  Absorption,  and  Rever- 
beration Noise 

The  complex  cepstrum  beamform  is  designed  to  minimize  the 
variance  of  log  normally  distributed  random  fluctuations  of  the  signal  due  to 
random  propagation  effects  and  to  some  extent  due  to  the  additive  log-normal 
noise  factor.  If  C.(f)  is  a random  propagation  effect  of  the  jth  channel  and 
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N.(f)  is  an  additive  zero  mean 
envelope,  then  the  jth  channel 


random  variate  with  a log  normally  modulated 
with  common  signal  spectrum  S(f)  is; 


where 


X.iO  = C.(f)  S(f)  + N.(f) 

.1  J J 

= S(f)  C.(f)  L.(f) 

J J 


(II-6) 


L.(f)  = 1 + C.(f)"^  N.(f). 

.1  J J 

For  visible  signals  it  is  expected  that  | C ,(f)|  ^ | N.(f)  | . To  the  extent  that  L.(f) 
can  be  approximated  as  a log  normal  variate,  beamforming  over  the  spatial 
distribution  of  channels  minimizes  transmission  effects  of  C.  and  pseudo- 
transmission effects  of  L..  The  variance  of  the  expected  log  amplitude  of  the 
signal  and  the  expected  phase  of  the  signal  due  to  such  independently  random 
transmission  effects  is  reduced  by  1 / J where!  channels  compose  the  complex 
cepstrum  beamform  of  the  signal.  It  is  noted  that  this  type  of  beamforming 
does  not  directly  reduce  the  expected  means  of  the  ambient  noise  as  does  time 
domain  beamforming.  But  indirectly  as  the  variance  of  the  log  normal  distri- 
bution is  lowered,  the  arithmetic  mean  amplitude  of  the  noise  will  also  be 
somewhat  reduced,  especially  the  occasional  large  amplitude  fluctuations. 


Before  beamforming  the  log  amplitude  of  data  in  the  noise 
buffer,  the  power  spectrum  of  noise  was  compared  to  that  of  the  signal.  Based 
on  stationarity,  it  was  expected  that  the  power  in  the  signal  buffer  should  ex- 
ceed the  power  in  the  noise  buffer  at  all  frequencies.  In  some  cases  due  to 
statistical  estimation  errors,  this  was  not  the  case.  In  those  cases,  the  log 
amplitude  of  the  signal  buffer  replaced  that  of  the  noise  buffer.  Thus,  the 
average  log  amplitude  of  the  data  in  the  noise  buffer  was  guaranteed  to  be 
equal  or  less  than  that  in  the  signal  buffer. 


To  correct  for  positive  bias  of  the  measured  mean  amplitude 
of  signal,  the  power  of  the  average  signal  was  reduced  by  the  power  of  the 
average  noise  as  follows,  where  the  subscripts  S+ N and  N refer  to  the 
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signal  and  noise  buffers,  respectively;  the  bar  refers  to  numerical  averaging 
over  channels:  and  the  E operator  refers  to  the  expected  value,  in  this  case, 
of  the  signal  amplitude  spectrum  L,og  A gl 

E(Log  Ag)  = Log  |Exp(Log  Ag^  ^)  - Exp(Log  Aj^)  ) . (II-7) 

This  expected  value  of  Log  Ag,  obtained  by  applying  a noise  bias  correction, 
must  be  positive  definite.  This  resulted  from  limiting  the  noise  power  esti- 
mate of  each  data  channel  to  power  less  than  or  ec;ual  to  the  power  of  data  in 
the  signal  buffer.  Tlie  real  part  of  the  signal  log  spectrum  was  replaced  by 
EtLog  -Ag)  • The  imaginary  part  of  the  signal  log  spectrum  remained  unchanged 
because  no  phase  bias  correction  can  be  derived  from  noise  power  information. 

At  a seismic  station,  the  observation  of  past  events  from  most 
regions  could  be  characterized  by  a nominal  v’alue  for  absorption,  T/Q. 

Events  from  most  regions  observed  at  NORSAR  were  observed  to  be  approxi- 
mately T/Q  = 0,  25,  as  shown  for  a number  of  events  by  flat  displacement 

-2  -3 

spectrum  near  the  corner  frequency  and  roll-offs  of  k or  k above  the 
corner  frequency.  It  is  expected  that  events  from  certain  localized  regions 
associated  with  high  heat  flow  could  be  observed  to  have  much  higher  absorp- 
tion factors,  T/Q.  Regional  absorption  corrections  were  applied  to  the  events 
following  the  method  of  Strick  (1970).  The  absorption-dispersion  operator 
used  to  correct  the  log  spectrum  is  as  follows: 

Uncorrected  for  absorption: 

Log  Z (kAf)  = Log  A (kAf)  + i (kAf)  . (II-8) 

S O O 

Corrected  for  absorption: 

Log  A=;=  (kAf)  = Log  A (kAf)  ‘ a-kAf(T/Q) 

= f/.^(kAf)  + 2kAf  (T/Q)  log(-i^^  ). 


f/)*g(kAf) 


(II-9) 


Before  applying  the  absorption  correction,  the  beamformed 
first  difference  of  the  phase,  A0g(kAf),  must  be  numerically  integrated  to 
I compute  the  phase.  The  initial  phase  at  k = 0 was  taken  as  zero.  The 

j ^ AOg(kAf)  was  exponentially  tapered  by  w"^  to  statistically  weight  the  influ- 

I ence  of  low  frequency  phase  information  and  to  minimize  roundoff  errors  due 

i 

to  the  integration.  The  taper  of  0.  96  used  was  determined  experimentally  by 
inspection  of  time  traces.  No  significant  distortion  of  the  time  trace  was 
noted  for  the  phase  angle  taper  greater  than  0.96.  The  integration  process 
■ used  to  compute  0 (kAf)  follows  as: 

Oj,(kA0  = w0g|(k-l)Af|  + A0g(kAf) 

k = r ‘ • K 0g(O)  = 0.  (II-IO) 

At  this  point  the  preceding  absorption  corrections  arc  applied  fo  the  real  and 
imaginary  part  of  the  log  spectrum.  The  absorption -di sper sio n corrected 
log  spectrum  of  the  signal  is  exponentiated, 

S(kAf)  = Exp  I Log  A":^(kAf)  + i 0g(kAf)j 

~ ^ <J)g(kAf)  -f  i sin  0g(kAf)]  , (II-ll) 

This  exponentiated  spectrum  is  transformed  to  the  time  domain; 

d>:<(nT)  = FFT'^  (s(kAO)  . 

to  produce  the  complex  cepstrum  beamformed  seismogram. 

As  a result  of  the  system  response  correction  applied  to  each 
channel  fed  into  the  above  complex  boamforming  process,  d*(nAT)  represents 
the  first  derivative  of  the  acceleration  of  earth  ground  motion. 

The  method  for  correcting  time  domain  data  for  system  re- 
sponse was  to  experimentally  fit  a physically  realizable  and  invertible  spec- 
tral operator  such  that  the  power  of  the  operator  matched  the  power  of  the 
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system  response  to  an  impulsive  ground  displacement  input.  The  spectrum 
of  the  seismic  system  response,  R(  Z) , is  given  as  follows,  where  Z ^ = 

exp(  - i 2 T7  kAf  T)  represents  the  spectrum  of  a time  delay  of  one  point  sampled 
every  T seconds: 


R(Z) 


Rj(Z) 

(poles) 


R^CZ) 

( zero's) 


( ^ 1 

3 1 

R(z)  = j n — ^ — - 
(m  "-“,=’■■>1 

1 

N 

i 

II 

(11-12) 


The  values  of  X and  were  determined  by  visual  curve  fitting  where  |R(Z)| 
fitted  the  measured  system  response  to  within  about  10  percent  in  the  range 
from  0.0  Hz  to  5,0  Hz,  Above  5.0  Hz  an  analogue  Butterworth  filter  was 
cascaded  into  the  system  response.  Since  we  could  not  obtain  consistent  spec- 
ifications of  the  roll-off  between  5.0  and  10.0  Hz,  no  attempt  was  made  to 
correct  the  response  in  that  frequency  range.  Note  that  R(  Z)  is  split  into  two 
parts;  the  first  part,  R^^(Z),  is  the  weighted  cascade  of  3 poles  with  roots 

Z = a.  and  weights  on  the  input  y..  The  second  part,  R ( Z) , is  the  cube  of  a 
J J ^ 

numerator  zero  factor  with  root  Z = 1.  The  inverse  of  R^(Z)  is  stable  and 

generally  increases  the  S/N  ratio  of  the  seismic  event.  This  operator  was 

applied  to  each  input  channel  was  follows,  where  d(nT)  is  the  input  with 

seismic  response  R(  Z) , and  d'(nT)  is  the  output  with  the  response  R^(Z): 


d'(nT) 


[d  (nT) 


Q.  d 
J 


((n  - 1)t)]/ V.  . 
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Forj=  1,  = 1.08,  and  = 0.7778,  The  output  of  this  first  stage  corre- 

sponding to  j = 1 replaces  the  input  d(nT)  to  form  a digital  cascaded  filter. 

The  two  following  operations  are  j = 2:  = 1.  39,  = 0.  0243;  and  j = 3: 

y = 3.  20,  Q = 0.  After  inverting  the  poles,  the  remainder  system  response 
/-!  3 , 

= (1  - Z ) is  a first  difference  operator  cascaded  three  times  or  a 


third  difference  operator.  Thus,  approximates  digitally  the  third 

derivative  of  ground  motion  displacement  or  first  derivative  of  ground  motion 
acceleration.  Since  = exp(-i  2 tr  kAf  T) , when  k = 0 (zero  frequency), 

the  magnification  of  (1  - Z”S'^  is  infinite.  For  this  reason,  the  inverse  is 
taken  approximately  to  avoid  blowing  up  low  frequency  noise  as: 


For  acceleration:  A(  Z) 


For  velocity:  V(Z) 


For  displacement:  D(  Z) 
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a Z 


-1 


u 
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aZ 
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aZ 


-1 
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1 - a z 

1 - v^Z 


-1 


-1 


1 - aZ 


-1 


1 - w^Z 


-1 


1 - aZ”^ 

^ ^ — (11-14) 
1 - w ^ Z 


The  numbers  used  in  the  above  spectral  operators  vere  a = 0,98,  u^  - 0.9, 

V = 0.91,  V = 0.90,  w = 0.92,  w = 0.01,  and  w = 0.90. 

1 Z 1 ^ -5 


The  time  domain  inverse  operators  are  applied  to  cem])ute  the 
acceleration,  velocity,  and  displacement  time  traces.  The  time  domain 
inverse  of  the  zeros  of  the  system  response  are  cascaded  by  the  denominator 
of  the  above  expressions  to  produce  acceleration,  velocity,  and  displacement 
wavelets,  valid  between  0.  25  and  5.0  Hz.  Inbetween  each  stage,  a shortpass 
cepstrum  filter  is  applied  to  improve  the  statistical  estimation  of  the  seismic 
ground  motion  and  to  remove  time  terms  longer  than  3,  2 seconds. 


The  smoothing  operator,  which  is  also  a short  pass  cepstrum 
operator,  was  performed  by  forming  a running  average  of  the  log  spectrum. 
The  method  used  to  average  was  to  use  staged  smoothing  with  binomial  coef- 
ficients. The  smoothed  log  amplitude  follows  as: 
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LogA"(kAf)  = 0.25 


2 Log  Ag  (kAf)  + Log  A^  |(k-j)  Afj 


+ Log  Ag  |(k+j)Af| 


(f)g(kAf)  = 0.25 


2 <f)  g (kAf)  + f/>g  |(k-j)Af  j 


(11-15) 


where  the  prime  indicates  the  state  of  the  data  before  smoothing;  unprimed  indi- 
cates the  state  after  one  stage  of  smoothing.  By  cascading  the  above  operator 
and  skipping  [0,1,2...]  points  after  each  stage  of  smoothing,  a Bartlett  window 
is  generated  to  smooth  the  log- spectrum.  For  example,  after  1 stage  the 
window  weights  are  ( 1 , 2,  1)  / 4;  after  2 stages  ( 1 , 2,  3,  2,  1)  / 9;  etc.  The  taper  of 
the  complex  cepstrum  (short  pass  operator)  for  the  Bartlett  window  produced 
by  I stages  of  smoothing  is; 

fj  cos^  i f i (.1-16) 

i=  1 


where  n is  the  number  of  point  delay  of  the  complex  cepstrum  and  N is  the 
number  of  data  points.  It  is  desirable  that  the  cepstrum  taper  be  relatively 
flat  up  to  at  least  32  points  delay  (sampled  at  0.05  seconds).  To  insure  this 
condition,  the  output  data  was  exponentially  tapered;  the  velocity  ground  mo- 
tion, exp  (n/0.99);  displacement  ground  motion,  exp(n/0.98).  The  resultant 

-M 

attenuation  of  an  M point  delay  operator  (Z  * ) by  smoothing  and  compen- 
sating with  the  exponential  taper  is  shown  in  the  following  Table  11-2. 

The  wavelets  at  delays  less  than  the  32-point  delay  shown  in 
the  table  are  progressively  closer  to  one  as  zero  point  delay  is  approached. 
It  appears  from  these  values  that  the  displacement  w.ivelet  echoes  can  be 
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sized  with  reasonable  accuracy  up  to  about  4,0  seconds.  This  should  be  ade- 
quate for  interpreting  pP  delays  of  explosions.  For  much  longer  delays,  for 
example,  due  to  multiple  explosions,  it  may  be  necessary  to  use  the  acceler- 
ation or  the  velocity  ground  motion  wavelet  instead  of  tlie  displacement  wave- 
let. The  Bartlett  window  used  to  smooth  the  above  wavelets  was  determined 
experimentally  as  a requirement  for  containing  the  dominance  of  long-period 
peaks  in  the  data,  ()ossibly  due  to  noise.  The  effect  of  the  smoothing  is  to 
statistically  average  out  such  long-period  peaks  over  a broader  band  of  fre- 
quencies. 


2,  Interactive  Data  Processing 


a.  Data  Conditioning 


In  this  phase  of  data  processing  (sec  Figure  II-8),  a seismic 
analyst  scans  a list  of  events  and  selects  one  for  tliscrimination  analysis.  The 
wavelet  is  automatically  corrected  for  absorption  based  on  the  parameter  T/Q, 
The  start  time  of  the  wavelet  is  picked.  Offset  of  the  initial  point  is  automati- 
cally corrected.  A taper  is  applied  to  avoid  hard  truncation  of  the  data.  The 
analyst  is  guided  bv'  displays  of  the  wavelet  and  the  real  cepstrum  of  the  wavelet. 

The  ABSORB  algorithm  corrects  for  values  of  T/Q  using  the 
prior  described  Strick  model.  The  LOGFFT  algorithm  is  as  follows; 


Xj  contains  128  points  of  data. 


contains  128  points  of  log  spectrum  in  multiplexed  FFT  for- 
mat with  imaginary  part  set  equal  to  zero. 


Log  |fFT(X^)^-^X^,  Rp{):  real  part  of, 

FFT’^  )V’^2^K^1*  display  (X^,Y^). 

After  examining  the  display  for  the  desired  start  point  of  the  signal,  the 
analyst  pushes  the  start  time  button  and  enters  the  number  of  points  delay  to 
start  the  wavelet.  Any  initial  offset  due  to  noise  is  automatically  removed. 
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FIGURE  II-8 
DATA  CONDITIONING 
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X(I)  = X(I)  - X(l)  * (N  - I + 1)/(N  - 1)  . 


(11-17) 


N is  automatically  taken  as  64  with  all  subsequent  operations  being  performed 
on  64  data  points.  If  desired,  the  analyst  may  taper  the  data  to  remove  abrupt 
truncations  of  the  signal's  coda.  A record  is  kept  of  the  number  of  times  that 
the  data  is  exponentially  tapered.  Later,  provision  is  made  to  invert  the  taper 
and  use  cepstrum  short  pass  filtering  to  remove  unwanted  later  arrivals  of 
energy.  Upon  exiting,  the  data  are  extended  with  zero's  out  to  128  points;  and 
the  real  and  imaginary  part  of  the  log-spectrum  is  saved  in  X^.  I’rom  this 
point  on,  X^  will  be  stored  on  disc  and  used  in  passing  the  wavelet  through  to 
the  next  subroutine.  The  replace  button  is  used  to  recover  an  erroneous  start 
time  or  taper  within  the  data  conditioning  subroutine. 

b.  Separate  Echo 

In  this  phase  (see  Figure  II-9),  the  selected  data  is  decon- 
volved to  remove  the  apparent  echo  which  the  analyst  detects  by  inspection  of 
the  real  cepstrum  display.  This  is  done  by  trial-and-error  input  of  a reflec- 
tion coefficient,  scattering  coefficient,  and  the  time  delay  in  data  points.  By 
pressing  a load  button,  the  echo  is  held  while  the  first  arriving  signal  log- 
spectrum  is  smoothed  to  short  pass  the  signal.  A second  press  of  the  load 
button  passes  through  the  echo  which  can  be  similarly  smoothed.  A third 
pass  of  the  load  button  reconstructs  the  signal  by  adding  the  first  arriving 
signal  to  the  echo  for  visual  comparison  with  the  original  data. 

The  algorithm  IFTEXP  displays  the  wavelet  and  the  real  cep- 
strum. 

fft'^  |Exp(X2)} -x^ 

FFT'^  {R^(X2)} 

DISPLAY  (X,  , Y ,) 

L 4 
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FIGURE  II-9 
SEPARATE  ECHO 


The  first  button  pushed  is  DECONVOLVE  which  follows-.  Let  y.  be  the  first 
arrivtil  signal  and  x the  input  data: 


''i  ' “l'’i-M-l  * " ’‘i 


aj  = RC  * SCl  a = RC  (I  - 2*SC)  . 


SC  is  between  0.  0 and  0,  5.  Given  0.  0 an  impulsive  echo  is  removed  by  the 
algorithm.  Given  other  values,  a symmetrical  multiple  echo  is  removed. 
Such  an  effect  would  shape  the  cepstrum  trough,  e.  g.  , a double  echo  or 
smoothed  echo.  The  analyst  can  try  any  desired  combination  of  delay  time, 
M;  reflection  coefficient,  RC;  and  scattering  coefficient,  SC,  until  the  cep- 
strum trough  is  apparently  removed. 

The  load  algorithm  stacks  the  data  and  echo  while  the  analyst 
smoothes  the  first  arrival  signal;  and  similarly  smoothes  the  echo.  After 
leaving  the  deconvolution  algorithm,  the  flag  IK  is  set  at  -1.  If  IK  >-l,  the 
deconvolve  algorithm  is  an  automatic  return  to  avoid  analyst  error.  After 
exiting  DECONVOLVE  the  data  are  stored  on  Z4,  and  the  deconvolved  signal 
is  stored  on  X , and  its  log-spectrum  on  X^-  Then  LOAD  is  entered: 

LOAD;  If  IK>2,  return 
IL  = 0 

IK  = IK  - I 
If  IK  ^ Z branch  to  (4) 

f X ^^1  (Reconstruct  data  by  adding  the 

branch  to  (2).  shortpasscd  signal  and  echo.) 


(4)  If  IK  = I branch  to  (I) 


7.  - X 7.  fecho  is  data  minus  sienall 


r 


(Data  in  Z^) 


IL  = IL  + 1 


IF  (IL.  EQ.  1)  Go  to  3 (return)  . 

This  LOAD  algorithm  automatically  saves  the  data,  deconvolved  signal,  and 
echo  as  they  are  generated  and  orders  the  data  as  follows:  After  first  press 

of  the  LOAD  button; 

IK  = 0 Z^:  original  data  minus  deconvolved  signal  (raw  echo) 

Z^'-  log  spectrum  of 

Z^:  original  data 

XjT  deconvolved  signal 

X log-spectrum  of  deconvolved  signal 

DISPLAY  (Z,,  X,  , Z,) 

4 1 1 

IK  = 1 X^:  log-spectrum  of  raw  echo 

Zy  original  data 

Z : smoothed  deconvolved  signal 

4 

Xjt  echo 

DISPLAY  (Zy  Zy  Xj) 


IK  = 2 Zy  original  data 

Zy  deconvolved  signal 
Zy  echo 

Z : deconvolved  signal  and  echo 

4 

Xj,;  original  data 
DISPLAY  (Z^,  Zy  Zy  Z^). 
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When  IK  = 0,  the  deconvolved  signal  is  smoothed  to  reduce  reverberation 


errors  and  short  pass  the  signal  estimate.  When  IK  = 1,  the  echo  estimate  is 
similarly  smoothed.  When  IK  = 2,  the  smoothed  signal  and  echo  are  added  to 
reconstruct  the  first  phase  and  its  echo.  Also,  the  original  dataare  stored  on 
X j to  prepare  input  to  the  next  subroutine  which  computes  the  residual. 

c.  Computation  of  the  Residual  Seismogram  to  Detect 
Multiple  Events 

The  signal  and  echo  of  the  first  phase  are  combined  to  model 
the  first  arrival.  This  is  subtracted  from  the  original  data.  Any  residual 
later  phase  left  on  the  record  can  be  picked  by  the  analyst  for  possible  subse- 
quent analysis.  The  ratio  of  the  removed  first  arrival  signal  to  the  residual 
is  calculated  from  zero  up  to  the  start  time  picked  for  the  possible  later  phase. 
This  S/N  ratio  is  displayed  along  with  the  possible  later  phase  and  its  cep- 
strum.  Also,  here,  the  correlation  between  the  signal  and  echo  of  the  first 
phase  is  computed.  If  desired,  the  later  phase  can  be  stored  on  the  disc  for 
the  next  entry  into  SPEED  and  thus  be  analyzed  as  a later  phase  of  the  same 
event. 

d.  Invert  Exponential  Taper  from  the  Signal  Estimate 

The  deconvolved  signal  estimate  is  corrected  for  the  exponential 
taper  applied  to  the  data  in  the  data  conditioning  described  in  subsection  a. 

This  involves  pressing  the  inverse  exponential  taper  button  until  the  number  of 
tapers  applied  to  the  data  arc  worked  off.  If  roundoff  errors  blow  up  the  data 
at  long  time  delays,  the  analyst  presses  the  zero  button,  inputting  the  number 
of  points  delay,  after  which  the  dataare  zeroed  out  to  eliminate  the  largest 
roundoff  errors.  During  the  sequential  inversion  of  the  exponential  taper, 
apparent  coda  noise  may  emerge  in  the  tail  of  the  signal.  This  can  be  removed 
by  pressing  the  smooth  button.  A replace  button  is  also  available  to  recover 
from  any  single  excessive  push  of  any  of  the  preceding  buttons.  After  each 
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press  of  a button,  the  current  state  of  the  deconvolved  signal  and  its  cepstrum 
are  displayed  to  control  the  next  button  push. 

e.  Absorption  and  Variable  Frequency  Magni tude 
Measurement 

The  variable  frequency  magnitude  method  of  discriminating 
earthquakes  and  explosions  is  based  on  magnitude  measurement  at  low  fre- 
quency and  high  frequency.  Explosions  should  be  characterized  by  signifi- 
cantly lower  magnitude  at  low  frequency  and  higher  magnitude  at  high  frequency. 
First,  the  analyst  is  asked  to  select  a low  frequency,  mid-band  of  frequencies, 
and  high  frequencies.  After  selecting  or  approving  the  pre-selected  values, 
the  log-amplitude  versus  frequency  spectrum  is  displayed.  An  option  is  pro- 
vided for  picking  any  two  frequency  points  and  linearly  extrapolating  the 
spectrum  either  to  the  selected  high  or  low'  frequency.  This  option  is  used  on 
events  only  if  the  selected  high  or  low  frequency  magnitude  appears  to  be 
dominated  by  noise.  A compute  magnitude  button  returns  a table  showing  the 
average  mid-band  log-amplitude,  the  difference  between  the  low  frequency 
magnitude  and  mid-band  magnitude,  and  the  difference  between  the  high  fre- 
quency magnitude  and  mid-band  magnitude.  Since  the  mid-band  frequencies 
were  selected  to  cover  the  expected  peak  frequency  of  seismometer  records, 
the  low  and  high  frequency  magnitude  measurements  can  be  approximated  as 
the  network  averaged  magnitude  of  the  event  as  seen  through  the  system  re- 
sponse plus  the  difference  magnitude  at  the  high  and  low  frequency,  respec- 
tively. 

Another  processing  option  can  be  performed  on  the  semi-log 
spvctrum  display  by  pressing  a button  for  determining  an  absorption  correc- 
tion, T/Q.  The  amplitude  spectrum  of  a signal  with  absorption  is  as  follows: 

Y(kAf)  = X(kAf)  Exp  (-1 /2  (T/Q)kAf)  (11-19) 
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where  T is  the  travel  time,  Q is  the  absorption  coefficient,  and  T/Q  is  ap- 
proximately constant  ♦’or  most  events.  The  parameter  T/Q  is  expected  to  be 
higher  than  normal  for  certain  high  absorption  paths.  Solving  the  above  equ- 
ation for  the  semi-log  spectrum,  LogX(kAf), 

Log  X(kAf)  = Log  Y(kAf)  i l/Z(T/Q)kAf.  (11-20) 

Rather  than  inputting  T/Q, which  is  unknown,  the  analyst  specifics  wlicther 
T/Q  is  positive  or  negative,  i'e  then  automatically  receives  a correction  for 
T/Q  ~ S,  wliich  produces  a small  rotation  of  the  spectrum  which  is  displayed 
for  the  analyst.  After  K pushes  of  the  T/Q  button,  the  analyst  may  obser''e 
that  part  of  the  spectrum  is  nearly  flat  and  therefore  the  absorption  correction 
is  completed.  At  this  point,  the  analyst  exits  or  selects  another  option  such 
as  measure  magnitudes,  and  the  value  of  the  T/Q  correction  is  held  as  anno- 
tation. Entry  back  into  condition  data  automatically  corrects  the  data  for 
absorption,  and  reduces  the  annotated  T/Q  back  to  zero.  Multiple  entries 
into  condition  data  are  taken  into  account  in  this  way.  In  addition,  the  cumu- 
lative T/Q  correction  of  the  basic  data  is  also  accounted  for  in  the  annotation. 

f.  Corner  Frcciuency  Measurement 

The  semi-log  spectrum  display  of  the  preceding  module,  cor- 
rected for  absorption,  is  converted  to  a log -amplitude  versus  log-frequency 
display  in  this  module.  By  picking  by  means  of  a cursor  two  frequencies  to 
the  left  of  a corner  frequency  and  two  frequencies  to  the  right  of  a corner  fre- 
quency, the  intersection  of  the  straight  lines  corresponding  with  each  pair  of 
points  yields  the  corner  frequency,  log-amplitude  of  the  corner  frequency,  and 
roll-off  to  the  left  and  right  of  the  corner  frequencies.  Up  to  four  corner  fre- 
quencies per  phase  can  be  computed  in  this  way.  Upon  exiting  this  module,  a 
table  of  corner  frequencies,  log-amplitudes,  and  roll-offs  are  displayed. 


F.  EXAMPLES  OF  HARDCOPY  FROM  SPEED 

The  Interactive  Seismic  Processing  System  (ISPS)  provides  an 
option  to  hardcopy  intermediate  results  of  processing.  For  the  case  of  the 
SPEED  module,  some  selected  hardcopy  is  described  and  shown  in  Figure 
II-IO,  pages  1 through  11. 

On  page  1 of  Figure  II-IO,  the  first  time  trace  is  a presumed 
explosion  from  Kazakh  corrected  for  instrument  response  to  approximate- 
ground  displacement  from  0.  3 to  5.0  Hz.  The  second  time  trace  is  the  posi- 
tive half  real  cepstrum  obtained  by  transforming  the  log-amplitude  spectrum 
into  the  time  domain.  Each  hashmark  on  the  cepstrum  trace  represents  a 
delay  of  four  points.  A possible  echo  appears  to  be  indicated  on  the  seventh 
point.  The  top  annotation  indicates  the  event  number,  source  location  and 
origin  time,  magnitude,  depth,  station,  channel  and  component  number,  sta- 
• tion  coordinates,  arrival  time,  event  azim-uth  and  distance,  and  sample  period. 

The  second  annotation  line  indicates  that  the  event  was  exponentially  tapered 
I twice,  started  on  the  first  point  of  input  data,  and  the  event  amplitude.  The 

last  line  c annotation  labels  the  real  cepstrum  of  the  first  phase  picked  and 
indicates  the  maximum  imiilitude  of  the  trace. 

On  page  2.  of  Figure  II-IO,  the  first  tra  v shows  the  displace- 
ment wavelet,  the  second  trace  shows  the  displacement  wavelet  with  the  pP 
echo  apparently  removed,  and  the  third  trace  shows  the  echo.  The  peak  am- 
I plitude  of  the  deconvolved  displacement  wavelet  is  increased  as  expected,  and 

one-sided  as  expected. 

On  page  3 of  Figure  II-IO,  the  deconvolved  wavelet  and  its 


cepstrum  are  shown.  The  second  line  of  annotation  shows  the  reflection  coef- 
ficients of  the  pP  echo,  the  scattering  coefficient  used  to  shape  the  i-cho,  num- 
ber of  points  delay,  and  amplitude. 


-t 


On  page  •},  the  deconvolved  displacement  wavelet  was  smoothed 
in  the  log-spectrum  domain  with  a Bartlett  window  of  five  points  (triangular 
half  -width).  This  short  passes  the  displacement  wavelet  as  shown  on  trace  I, 
The  cepstrum  of  the  short  passed  displacement  wavelet  is  shown  on^trace  2. 

On  page  5,  the  first  time  trace  shows  the  displacement  wavelet. 
The  second  time  trace,  the  short  passed  deconvolved  displacement  wavelet. 
The  third  trace,  the  short  passed  echo.  The  fourth,  the  displacement  data 
reconstructed  by  adding  the  short  passed  deconvolved  wavelet  to  the  short 
passed  echo.  Comparison  of  the  reconstructed  data  and  original  data  shows 
close  correspondence.  The  amplitude  of  the  reconstructed  wavelet  is  down 
50  nm  or  about  1 'Ti  from  the  original  data. 

On  page  6 of  Figure  II- 10,  the  difference  between  the  original 
data  and  reconstructed  data  is  sliown  on  trace  1.  The  maximvim  difference  is 
658  nm  or  about  12%  of  the  maximum  signal  am})litude.  Thc'  largest  errors 
occur  before  the  start  of  the  signal,  under  the  echo,  and  at  the  end  of  thc  data 
trace.  Under  the  first  arrival  pulse,  the  error  is  much  less,  several  percent, 
Thc  second  trace  shows  the  cepstrum  of  thc  residual  trace.  The  last  line  of 
annotation  shows  the  correlation  coefficient  between  the  signal  and  echo  as 
-0. 988. 

On  page  7,  the  rt.-sult  of  using  a cursor  to  pick  a possible  later 
phase  is  sh  >un.  The  first  trace  shows  thc  later  phase,  the  second,  its  cep- 
strum. The  second  line  of  annotation  indicates  that  the  later  phase  was  time 
delayed  11  points  and  its  maximum  amplitude  is  514  nm.  The  bottom  line  of 
a otation  indicates  a signal  to  noise  ratio  of  12.4  db  based  on  power  of  thc 
data  compared  to  the  power  of  the  residual  trace  obtained  by  removal  of  thc 
reconstructed  data,  comjjuted  up  to  the  presumed  arrival  of  the  later  phase. 
The  so-called  later  phase  was  picked  for  demonstration  purpose  and  is  not 
consider<'fl  a v'alirl  later  phase. 
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On  page  8 of  Figure  II-IO,  the  short  passed  deconvolved  wave-  . 

let  is  shown  with  the  applied  exponential  taper  inverted.  The  amplitude  of  the  j 

signal  was  increased  from  4,  560  nm  to  6,420  nm.  This  is  the  best  measure-  f 

ment  of  signal  amplitude.  The  cepstrum  is  shown  on  the  second  trace.  Small 
positive  secondary  delays  are  indicated  with  4-point  and  7-point  delays.  These 
can  be  seen  on  trace  1,  the  displacement  first  arrival  wavelet,  as  possible 
small  secondary  arrivals  as  peaks  with  4-point  and  7-point  delays,  consistent 
with  the  cepstrum.  i 

On  page  9,  the  log -amplitude  spectrum  is  plotted  as  a function 
of  frequency.  The  vertical  lines  show  on  the  left  the  low  frequency  used  to 
compute  low  frequency  magnitude.  The  two  middle  lines  show  the  two  fre- 
quencies for  averaging  mid-band  magnitude.  The  right  vertical  line  shows 
the  frequency  used  to  compute  high  frequency  magnitude.  The  T/Q  factor  of 
zero  indicates  no  absorption  correction  was  applied  to  the  input  data.  The 
bottom  line  shows  smoothes  equal  to  zero,  indicating  that  the  log-amplitude 
was  not  smoothed.  The  log -amplitude  of  2.  92  - 4.  47  shows  the  range  of  log- 
amplitude  in  the  plot,  with  base  10  logarithms. 

On  page  10,  the  s erhi -logarithm  spectrum  of  the  first  arriving 
signal  pulse  is  shown  corrected  for  absorption  using  a T/Q  equal  to  0.  225. 

The  log -amplitude  range  has,  of  course,  been  narrowed  by  the  absorption 
corretion  from  4.00  to  4.  50.  The  absorption  correction  will,  as  expected, 
significantly  increase  the  magnitude  of  the  event.  The  average  log-amplitude 
between  1. 00  and  1.75  Hz  is  4.  37.  The  difference  log -amplitude  at  0.  3 Hz  is 
0.  13  magnitude  units.  The  difference  log-amplitude  at  2.  50  Hz  is  -0.03  mag- 
nitude units. 

On  page  11,  the  log -amplitude  spectrum  of  the  deconvolved 
displacement  wavelet  is  plotted  against  the  log -frequency  from  0.  45  Hz  to  3.  5 
Hz.  This  plot  includes  the  previous  correction  for  absorption.  The  integers 
I and  2 next  to  the  plot  indicate  where  corner  frequencies  were  picked.  The 

♦ 
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log-frequency  from  -0.  35  to  0.  55  indicates  the  log-frequency  range  of  the 
horizontal  axis, and  the  specification  amplitude  from  3.  98  to  4.  48  indicates 
the  range  of  log-amplitudes  on  the  vertical  axis.  The  last  four  lines  of  anno- 
tation provides  a table  of  corner  frequency  information.  The  corner  frequen- 
cies for  pick  1 and  pick  Z are  giv>^en  under  column  1 and  column  2,  respectively, 
in  log  base  ten  units.  The  magnitude  at  the  corner  frequencies  is  also  given 
in  log  base  ten  units.  The  left  slope  indicates  the  power  f^  roll-off  to  the  left 
of  the  corner  frequency  pick.  The  right  slope,  to  the  right  of  the  correspond- 
ing corner  frequency  peak.  The  results  for  this  event  indicated  a corner  fre- 
quency at  2.4  Hz  nearly  flat  to  the  left  of  the  corner  with  roll-off  of  1. 8 power 
to  the  right  of  the  corner.  Another  corner  was  picked  at  2.9  Hz  with  roll-off 
of  3.0  power  to  the  right  of  that  corner.  Sax  (1975)  indicated  that  measure- 
ments of  log-amplitude  and  corner  frequencies  could  be  used  to  discriminate 
between  earthquakes  and  explosions.  In  that  study,  maximum  entropy  spectral 
estimates  were  used.  In  this  studv,  conventional  discrete  Fourier  transforms 
were  linearly  interpolated  at  uniformly  coupled  log-frequencies  after  correct- 
ing for  absorption,  based  on  the  assumption  that  at  least  part  of  the  displace- 
ment wavelet  sp'"Ctrum  is  flat. 
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SECTION  III 


EVALUATION  OF  ROUTINE  DISCRIMINATION 
PROCESSING  PROCEDURES 

A.  INTRODUCTION 

In  the  context  of  a real  time  surveillance  system  environment, 
it  is  estimated  that  an  average  of  fifty  events  per  day  will  be  collected.  After 
each  event  is  located,  the  retrieved  station  waveforms  are  time-aligned  to 
measure  accurate  focal  parameters  and  to  measure  accurate  station  param- 
eters such  as  station  magnitude  and  arrival  time  deviations,  first  motion 
sense,  and  other  station  parameters.  The  station  waveforms  are  further 
reduced  to  one  or  several  event  waveform  representations  for  routine  dis- 
crimination processing.  One  of  the  possible  reduction  methods  used  is  cep- 
strum  beamforming,  which  reduces  ambient  and  coda  noise  subject  to  the 
constraint  that  log-amplitude  spectrum  of  the  estimates  source  signal  is  the 
average  of  all  the  input  channels  fed  into  the  cepstrum  beamforming  process. 
The  algorithm  used  to  compute  the  cepstrum  beamforming  data  is  described 
in  the  preceding  section. 

The  discrimination  processing  performed  on  data  collected  in 
a real  time  surveillance  mode  must  all  be  done  with  sufficient  speed  to  keep 
up  on  an  average  with  the  event  load  collected  by  the  system.  This  means 
that  all  routine  discrimination  processing  must  be  done  in  less  than  one-half 
hour.  This  provides  sufficient  event  processing  capacity  to  work  off  back- 
logs due  to  occasional  swarms  of  events  in  a reasonable  span  of  time. 

For  the  purpose  of  quickly  evaluating  proposed  standards  of 
short-period  discrimination  processing,  a data  base  of  thirty-five  events 
was  prepared.  The  data  base  consisted  of  twenty  earthquakes  and  fifteen 
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presumed  explosions.  All  events  were  selected  in  the  teleseismic  range 
from  25  degrees  to  85  degrees  epicentral  distance.  Two  classes  of  earth- 
quakes were  selected.  In  one  class  labelled  as  simple  earthquakes,  the 
events  were  impulsive  with  the  maximum  peak  occurring  in  the  first  several 
seconds.  In  the  other,  labelled  as  complex  earthquakes,  the  maximum  peak 
occurs  later  than  the  first  several  seconds.  Earthquake  events  were  as- 
signed to  these  groupings  to  equally  weight  earthquake  results  for  simple 
and  complex  events  and,  if  possible,  differentiate  results  obtained  for  these 
two  groupings.  The  magnitudes  of  events  selected  for  the  data  base  nearly 
uniformly  cover  the  range  from  4.  5 to  6.  0 for  all  classes  of  events:  simple 

earthquakes,  complex  earthquakes,  and  presumed  explosions.  All  earth- 
quakes selected  were  shallow  focus  where  determined  or  otherwdse  of  indi- 
cated normal  depth  (33  km).  The  presumed  explosion  data  included  twelve 
events  from  the  eastern  Kazakh  region,  including  two  presumed  FNE's;  two 
events  from  NTS;  and  one  from  Colorado.  All  the  data  used  were  recorded 
at  NORSAR  on  twenty-two  subarrays.  All  the  data  preparation  and  inter- 
active processing  were  performed  starting  at  the  level  of  1 32  individual 
sensor  records  and  working  dov^'n  to  the  level  of  the  several  event  waveform 
representations  used  for  discrimination  processing.  The  processing  proce- 
dures developed,  when  perfected,  could  be  adapted  to  p..>cessing  a sizeable 
network  of  sensors  within  the  time  constraints  imposed  by  the  requirement 
for  real  time  rate  event  classification  processing. 

The  purpose  of  the  data  base  is  to  provide  preliminary  screen- 
ing of  the  effectiveness  of  discrimination  methodologies  which  are  possible 
candidates  for  routine  on-line  discrimination  processing.  Although  the 
screening  procedure  should  quickly  eliminate  discrimination  procedures 
which  perform  poorly,  any  positive  results  should  be  further  tested  with  a 
much  larger  data  base  before  acceptance  of  the  procedure  is  finalized.  A 
description  of  the  earthquake  event  information  is  shown  in  Table  III-l;  of 
the  presumed  explosion  data,  in  Table  III-2. 
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TABLE  III- 


'■  Presumed  PNE. 


On  Figure  UI-1,  plots  are  shown  of  the  representative  event 
waveforms  obtained  by  cepstrum  beamforming.  The  waveforms  are  shown  as 
seen  through  the  system  response  and  as  corrected  for  ground  motion  (accel- 
eration, velocity,  displacement  as  specified  by  the  user)  in  the  band  of  fre- 
quencies between  0.  25  and  5.0  Hz. 

B.  BASELINE  FOR  SHORT-PERIOD  DISCRIMINATION  PROCESSING 

In  Section  II,  an  example  of  the  application  of  short-period 
interactive  processing  was  shown.  In  the  example,  it  was  possible  with  very 
little  ambiguity  to  use  the  cepstrum  to  deconvolve  an  apparent  echo  to  derive 
the  expected  one-sided  displacement  wavelet  and  echo  with  correlation  coeffic- 
ients close  to  unity.  However,  in  attempting  to  repeat  the  performance  on  all 
the  events  in  the  screening  data  base,  it  was  no  longer  found  to  be  possible  to 
consistently  and  unambiguously  use  the  cepstrum  to  deconvolve  an  apparent 
highly  correlated  echo  as  in  the  example.  To  develop  a satisfactory  proce- 
dure for  routine  discrimination  processing,  a baseline  procedure  will  be  de- 
fined. The  procedure  will  be  run  through  the  screening  data  base  to  establish 
performance  figures  for  the  baseline  procedure.  By  systematically  evolving 
procedures  which  can  be  shown  to  perform  more  effectively  than  the  baseline 
procedure,  a satisfactory  procedure  for  discrimination  processing  can  even- 
tually be  developed.  The  interactive  SPEED  program  is  a useful  research 
tool  for  this  purpose.  By  specifying  a processing  procedure,  the  entire 
screening  data  base  can  be  run  in  two  or  three  hours  to  provide  basic  data  for 
preliminary  evaluation  of  the  procedure. 

The  procedure  followed  to  establiah  a baseline  for  discrimina- 
tion was  to  measure  variable  frequency  magnitude  data  on  seismic  data  cor- 
rected for  system  response  to  ground  displacement  and  ground  acceleration. 

All  the  event  data  were  processed  as  follows; 
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The  displacement  data  of  each  event  were  contained  in  a 6.  4 
second  time  window.  The  data  were  tapered  by  exp(-0.01  t)  to  increasingly 
weight  the  data  at  the  start  of  the  signal  and  to  reduce  discontinuities  at  the 
end  of  the  data  sample.  The  log  spectrum  of  the  data  was  smoothed  with  a 
Bartlett  window  (0.  1 1 , 0.  Z2,  0.33,  0.22,  0.11).  After  transforming  back  to 
the  time  domain,  the  applied  taper  was  inverted  by  applying  the  inverse  expo- 
nential taper  exp  (+0.01  t)  to  the  data.  The  low  frequency  magnitudes  were 
measured  at  0.  3,  0.  4,  and  0.  5 Hz.  The  high  frequency  magnitudes  were 
measured  at  2.5,  3.0,  and  4.0  Hz.  The  center  band  magnitudes  were  av’or- 
aged  from  0.7  5 Hz  to  1.7  5 Hz.  The  acceleration  wavelet  of  each  event  was 
contained  in  a 3,  2 second  time  window.  The  data  were  tapered  by  exp(-0.  02  t). 
The  same  Bartlett  window  was  applied.  The  same  frequencies  were  measured 
with  additional  high  frequencies  measured  at  5.0,  6.0,  and  8.0  Hz.  No  sub- 
jective measurements  of  echo  delays,  deconvolution,  or  absorption  correc- 
tions have  been  applied  to  the  data  for  this  baseline  procedure.  In  future  de- 
velopments, such  procedures  could  be  applied,  but  only  if  significant  gains 
over  the  baseline  results  can  be  expected. 


C.  RESULTS  OF  BASELINE  DISCRIMINATION  PROCESSING 
1.  Magnitude  from  Spectrum  Measurement 

The  average  log  amplitude  measurements  between  1.0  and  1.75 

Hz  are  related  to  network  magnitude  estimates,  m,  , as  follows: 

bs 

For  displacement  log  amplitudes: 


m,  = log  D + 1.  353  log  sinA  + 1.993 
bs 


For  acceleration  log  amplitudes: 
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m,  = log  A + 1.  353  log  sinA  + 2.77 
bs 
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where 


D - logarithm  of  the  displacement  spectral  amplitude 
averaged  between  1.00  and  1.7  5 Hz 

^ = epicentral  distance  in  degrees 

A = logarithm  of  the  acceleration  spectral  amplitude 
averaged  between  1.00  and  1,7  5 Hz. 

The  constants  1.  353,  1.993,  and  2.  77  were  obtained  from  the  data  by  least 
squares  regression  of  the  data.  The  results  are  shown  on  Figure  III-2.  For 
all  events,  the  standard  deviation  between  computed  and  world-wide  network 
magnitudes  was  0,  27  for  magnitudes  computed  from  average  log  displacement 
amplitude  and  0.  21  for  magnitudes  computed  from  average  log  acceleration 
amplitude.  In  the  latter  case,  three  anomalously  low  magnitudes  (about  one 
magnitude  unit)  were  computed  using  acceleration  data.  This  occurred  with 
three  earthquakes:  events  38,  55,  and  89.  The  result  was  probably  due  to 

nodes  in  the  radiation  pattern  causing  an  observed  hole  in  the  spectrum  at 
less  than  1.75  Hz.  The  computation  of  magnitude  from  displacement  ampli- 
tude data  for  these  events  was  also  low  but  by  less  than  0.  5 magnitude  units 
due  to  heavier  weighting  of  amplitudes  at  the  1.0  Hz  end  of  the  range,  which 
are  frequencies  lower  than  the  observed  hole  in  the  spectrum. 

The  lower  standard  deviation  of  rriagnitude  measurements  from 
log  acceleration  amplitude  measurements  is  probably  due  to  the  fact  that  the 
maximum  peak-to-peak  amplitude  through  the  system  response  corresponds 
more  closely  to  ground  acceleration  than  to  ground  displacement.  Based  on 
source  theory,  the  sealing  of  magnitude  as  the  square  root  of  . nergy  requires 
that  both  the  log-displacement  amplitude  and  corner  frequency  be  used  to 
estimate  magnitude  from  spectrum.  It  is  possible  that  the  standard  deviation 
of  magnitudes  computed  from  ground  displacement  can  be  further  reduced  by 
such  scaling. 
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che  preceding  section,  it  was  shown  that  network  magnitude 
of  the  events  could  be  reasonably  well  estimated  by  averaging  log  spectrum 
measurements  between  1.0  Hz  and  1,75  Hz.  In  order  to  statistically  compare 
low  frequency  and  high  frequency  magnitude  measurements  of  events  of  differ- 
ent magnitude,  the  average  magnitude  between  1,0  and  1.75  Hz  was  sub- 
tracted from  the  low  frequency  magnitude  and  high  frequency  magnitude 
measurements.  These  variable  frequency  magnitude  differences  are  shown 
in  Table  III -3. 

The  following  are  definitions  of  column  headings  used  in  Table 
III-3  and  other  tables: 


• °M 
1.  75  Hz. 

• M(0.  3)  - 

• M(0.  4)  - 

• M(0.  5)  - 

• M(2.  5)  - 

• M(3.  0)  - 

• M(4,  0)  - 

• y(0.3)  - 

• y(o.  4)  - 

• V(2.  5)  - 

• y(3.0)  - 

• Log  A.  . 


- average  log  displacement  amplitude  between  1.0  and 

log  displacement  amplitude  at  0.  3 Hz  minus  Log 
log  displacement  amplitude  at  0,  4 Hz  minus  Log 
log  displacement  amplitude  at  0,  5 Hz  minus  Log 
log  displacement  amplitude  at  2.  5 Hz  minus  Log 
log  displacement  amplitude  at  3.0  Hz  minus  Log 
log  displacement  amplitude  at  4.  0 Hz  minus  Log 
log  log  spectrum  rolloff  between  0.  3 and  0,  4 Hz. 

log  log  spectrum  rolloff  betv  een  0.  4 and  0.  5 Hz. 

log  log  spectrum  rolloff  bet\  een  2.  5 and  3.0  Hz. 

log  log  spectrum  rolloff  between  3.0  and  4.0  Hz, 

- average  log  acceleration  given  for  reference  only. 
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M = spectral  magnitude  levels  7 = log  log  spectral  slope 


1 


Positive  values  of  the  y's  indicate  increasing  amplitude  with  frequency;  nega- 
tive, decreasing.  Multiplication  of  the  y's  by  6.0  would  indicate  measured 
dB/octaves  on  a log  amplitude  versus  log  frequency  plot  of  the  spectrum. 

Thus  the  statistics  in  Table  III-3  characterize  the  high  and  low  frequency 
spectrum  by  levels  and  slopes  of  the  log  log  spectral  plot  of  the  event.  The 
procedure  used  to  generate  Table  III- 3 was  standardized.  Several  of  the 
events  were  run  several  times;  the  results  were  reproducible.  At  the  end  of 
Table  III-3,  the  MEAN  (Q)  represents  the  average  earthquake  spectrum  sta- 
tistics; MEAN  (X),  average  of  central  Asia  presumed  explosions  (excluding 
event  60  as  an  outlier);  and  MEAN  (N) , average  western  United  States  pre- 
sumed explosions  (excluding  event  59  as  an  outlier).  SD(Q)  is  the  estimate 
of  the  standard  deviation  of  the  representative  earthquake  population.  In  the 
western  United  States  presumed  explosion  population,  event  59  was  cast  out 
from  the  mean  as  an  outlier,  as  was  event  60,  from  the  central  Asian  popula- 
tion. 

A special  case  for  application  of  the  event  spectral  discrimina- 
tion measurements  is  to  plot  the  low  frequency  magnitude,  m , versus  the 

high  frequency  magnitude,  m , where 

H 

m.  = m_  + M(0.  3) 

and 

m = m + M(4. 0).  (III-3) 

H D 

The  results  are  shown  on  Figure  III-3. 

3.  Multivariate  Spectral  Discrimination 

The  variable  frequency  magnitude  data  (see  Table  III- 3)  are 
transformed  into  detectability  data.  The  measurement  vector; 


X ; [M(0.3),  M(0.4),  M(0.5),  M(2.5),  M(3.0),  M(4.0), 

y(0.3),  y(0.4),  y(2.  5),  y(3.0)]  (111-4) 


■f  - Central  Asia  Presumed  Explosions 
^ - Western  U.  S.  Presumed  Explosions 
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has  subtracted  from  it  the  mean  of  the  earthquake  population  Eq(X), 


whei 


Eq(2i) 


[0.79,  0.93,  0.  92,  -1.12,  -1.54,  -2.  02,  1.  20, 

-0.  10,  -4.  50,  -3.  12]  , (III-5) 


and  is  divided  by  the  standard  deviation  of  the  earthquake  population, 
where: 


SD^(X), 


sDq(X)  : 


[0.41,  0.35,  0.27,  0.  37,  0.41,  0.41,  0.61,  1.02, 

1.82,  2.31]  . (111-6) 


For  earthquakes,  the  detectability  is  a unit  normal  Z-statistic.  The  mean 
y(2.  5)  and  '>'(3.  0)  were  taken  at  the  50  percent  levels,  and  the  standard  devia- 
tions were  computed  by  weighting  positive  deviations  only  which  are  most  rel- 
evant for  discrimination.  Obtaining  X from  Table  III- 3 and  using  the  above 
values  for  E^fX)  and  SD^(X) , detectabilities  were  computed  for  all  35  events 
as  |x  - Eq(2£))/SD^(X)  with  the  results  shown  on  Table  II1-4,  for  the  earth- 
quake and  presumed  explosion  data. 


The  information  contained  in  Table  111-4  is  the  basic  data  for 
designing  an  optimum  least  squares  multivariate  discrimination  process  using 
measures  of  spectral  levels  and  spectral  rolloffs.  Each  component  of  the 
earthquake  detectability,  in  Table  1II-4  is  theoretically  a unit  variance, 

zero  mean,  Gaussian  statistic.  Each  component  of  one  or  more  presumed 
explosion  populations  is  possibly  detectable  as  a large  variation  from  the 
mean.  For  example,  the  mean  detectability  of  the  central  Asian  presumed 
explosion  population  is: 

: [-0.34,  -0.  42,  -0.57,  2,05,  2.04,  1.66,  -0.25, 

-0.  08,  0.  10,  -0.  84]  , (ni-7) 

The  dominant  statistics  are  significantly  elevated  high  frequency  spectral 
levels.  Low  frequency  levels  tend  to  be  consistently  lower  than  earthquakes, 
and  the  rolloff  between  3.  0 and  4.  0 Hz  is  steeper  than  earthquakes.  The 
mean  detectability  for  western  United  States  presumed  explosions  is: 
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spectral  magnitude  levels  Y - log  log  spectral  slope 


DETECTABILITY  DATA 
(PAGE  2 OF  2) 


F 


[ 0.34,  0.  89,  1.68,  -0.40,  0.30,  1.59,  2.06,  1.  45, 
1. 43,  1.  53] . (III-8) 


The  dominant  statistics  are  progressively  increasing  spectral 
levels  at  low  frequencies,  significantly  elevated  spectra  at  4.0  Hz,  signifi- 
cantly more  steeply  increasing  spectral  slopes  at  low  frequencies,  and  less 

rolloff  at  high  frequencies.  The  coherence  between  E(D  ) and  E(D  ) is 

"“X  N 

0.05  indicating  that  the  events  represented  by  the  average  detectability  E(D  ) 
are  detectable  to  an  almost  completely  independent  basis  from  those  events 
detectable  by  E(Dj^).  Thus  at  least  two  distinct  sets  of  weights  must  be  de- 
rived to  optimally  detect  presumed  explosions  from  western  United  States  and 
from  central  Asia. 


The  performance  of  an  optimum  m.ultivariate  detector  depends 
not  only  on  the  mean  of  the  signal  classes  being  detected,  but  also  on  their 
variation  from  the  mean  signal  model.  The  difference  from  mean  signal  mod- 
els E(D^)  and  E(Pj^j  are  shown  on  Table  III-5.  The  covariance  matrices 
for  earthquake  errors  and  presumed  explosion  errors  are  shown  on  Table 
III-6.  The  large  variances  for  y(2.  5)  and  y(3.  0)  (the  last  two  diagonal  values 
of  earthquake  area  variances  on  Table  III-6)  are  due  to  using  positive  devi- 
ations to  compute  the  last  two  elements  of  SDq(X). 


The  least  squares  equation  for  optimizing  the  detection  of 
presumed  explosions  from  central  Asia  is  as  follows.  Elements  of  the  sym- 
metric earthquake  correlation  matrix,  [ R]  , are  R ; 

ij 


where 


R. . 
ij 


1 

K - 1 


Z 


X.(k)  X.(k) 
1 J 


(III -9) 


K = the  number  of  earthquakes  used,  and 

X^  = the  ith  element  of  X,  which  is  the  vector  of  measured 
discriminants. 
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TABLE  III-6 

COVARIANCE  MATRICES 
(PAGE  1 OF  2) 

Earthquakes  Error  Data 
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TABLE  III-6 
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0.  08 


COVARIANCE  MATRICES 
(PAGE  2 OF  2) 


Presumed  Explosions  Error  Data 
Central  Asia  Events 
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The  summation  K is  over  the  earthquake  population.  The  optimization  is 
furthered  by  casting  out  large  negative  deviations  of  the  dot  product  ‘X  | 

<-  SQRT  (D^)  • E Thus  the  least  squares  process  makes  no  attempt 

to  reduce  the  variance  contributions  of  large  negative  earthquake  deviations  as 
these  are  easy  to  discriminate.  If  a much  larger  data  base  were  used  to  de- 
sign the  optimum  multivariate  discriminant,  then  |e  (D^)  ' X |<  0 would  be 
cast  out  as  only  positive  earthquake  errors  result  in  false  positive  earthquake 
identification.  As  a result  of  this  condition,  the  large  negative  detectabilities, 
which  obviously  indicate  earthquakes,  are  unconstrained. 

A covariance  matrix  [ R']  is  constructed  for  a combined  popu- 
lation of  earthquakes  and  presumed  explosions  from  central  Asia,  Let 

U (m)  be  the  vector  of  average  detectabilities  of  presumed  explosions  from 
the  mth  region  or  explosions  medium  type.  For  example,  from  the  central 
Asian  set,  m = 1,  and  U.(l)  is  the  ith  element  of  U(l),  elsewhere  referred 
to  as  E (D^); 

K 

[S]  = S..{m)  = ^ Y]  [U.(k,m)  - ] 

m ij  K - 1 * ^ 1 1 

k:=  1 

* [U,(k,  m)  - J.(m)]  . (III-IO) 

J 

[ S]  is  the  covariance  matrix  of  signal  errors  from  the  mth  region  or  explo- 
m 

sion  medium  type.  [S]  is  shown  or.  the  bottom  of  Table  III-6  for  m = 1, 

m 

central  Asian  presumed  explosion  types,  and  for  m = 2,  western  United  States 
presumed  explosion  types.  The  covariance  [ R']  which  combines  the  variance 
due  to  earthquake  errors  and  presumed  explosion  errors  is: 

R'..(m)  = a • [ U.(m)  U.(m)  + S.  .(m)  ] + { 1 -a)  R. . . (III-ll) 

U 1 J ij  U 
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The  scalar  factors,  a,  where  (0<all)  and  (1-a),  provide  the  designer  the 
option  to  relatively  weight  the  minimization  of  errors  due  to  presumed  explo- 
sion deviations  and  the  errors  due  to  earthquakes  in  the  assumed  mixed  pojni- 
lation.  For  example,  if  a=  1 , the  multivariate  discriminant  minimizes  the 
\'ariance  of  the  explosion  population  and  ignores  the  variance  <lue  to  earth- 
quakes. If  a=  0,  presumed  explosion  errors  are  neglected  and  positive  devi- 
atioi  s of  earthquakes  are  minimized.  Several  trials  have  indicated  that  the 
multivariate  discriminants  are  more  effectivt'  if  a is  nlose  to  zero. 

The  least  scpiares  solution  for  the  opiimuir  omltiva  riati'  dis- 
criminant W(m)  f<tr  presumed  explosion  events  from  the  mth  ri'gior  is: 

W(m)  = (III-12) 

where  the  elentenls  of  the  vector  l-I(m)  , the  expec  ted  '.-alue  of  detectabilities 
of  the  mth  signal  region  or  mediiim  type,  is  U,(m). 

In  obtaining  stable  solutions  for  W(l)  for  the  central  Asian 
events  and  W(2)  for  the  western  United  States  events,  the  off-diagonal  ele- 
ments of  the  matrix  [ R'j  must  be  nniltiplied  by  a positive  factor,  h,  equal  to 
■O'  less  than  uie.  The-  solution  is  eonsidired  stable-  if  the  maximum  ab^.olut■ 
change  in  the  solution  (given  an  iteration  of  tlie  off-diagonal  weighting  factor 
I))  is  less  than  t . (1-b)  • SORT  ^W(n)  • W(nlj.  In  our  case,  h was  taken 
as  0,9^1  and  r as  3.0.  The  solution  vectors  W(n)  wore  normalizeel  so  that 
the  expected  value  of  a presumed  explosion  was  l.O  and  of  an  t-arthquake  was 
0.0.  For  the  central  Asian  presumed  explosions,  the  off-diagonal  elements 
of  [R']  were  multiplied  by  0.99  and  for  the  western  United  States,  0.93.  This 
satisfied  the  above  conditions  for  a stable  solution.  The  solution  for  the 
central  Asian  presumed  explosions  is: 

W(l)  = [-0.07,  -0.07,  -0.07,  0,20,  0.10,  0.20,  0,05,  0.09, 
0.00,  0.09]  (III-13) 
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and  for  the  western  United  States  presumed  explosions, 


W(2)  = [0.02,  0.03,  0.13,  0.05,  0.02,  0.08,  0.  05,  0.26, 

0.  03,  0.  08]  . (III-14) 

The  coherence  between  W(l)  and  W{2)  is  0.31.  Thus,  W(l) 
and  W{2)  provide  a nearly  independent  basis  for  detecting  events  from  the 
two  regions.  Figure  III-4  shows  the  results  obtained  by  plotting  W(l)  • X, 
the  componv;it  directed  toward  the  central  Asian  set  of  presumed  explosions 
versus  W(2)  ■ X,  the  componen*'  directed  toward  the  western  United  States 
presumed  explosion  component. 

The  set  of  vectors  W(n),  which  are  used  to  weight  discrimi- 
nants of  events  from  N regions  or  medium  types,  are  in  general  correlated. 
To  obtain  an  independent  and  orthogonal  basis  for  detecting  events  of  N types, 
the  set  W(n)  is  replaced  by  a set  of  unit  normalized  orthogonal  vectors,  W'(n). 
Let: 

Y(n)  = W(n)/|W(n)l,  n=l,2,...N.  (III-15) 

The  correlation  matrix  of  the  set  of  vectors  Y(n)  is  an  N x N matrix  [ A]  , 
wher  e: 

= Y(k)  • Y(n)  . (III-16) 

From  the  above  nor mali  zation,  the  diagonal  elements  of  [a]  (i  = j)  arc  1.0; 
the  off-diagonal  element  magnitudes,  the  coherence  between  the  vectors  Y (k) 
and  Y(n);  and  the  sign  of  the  off-diagonal  elements,  the  sense  of  the  correla- 
tion between  Y(k)  and  Y(n).  The  solution  for  the  matrix  of  independent  orthog- 
onal unit  vectors  W'(n)  are  as  follows: 

[ W'(n)  ] = [ A]  W(n)  ] , (III-17) 


111-29 


Central  Asia  Presumed  Explosion  Region 


•V 


MULTIVARIATE  DISCRIMINANTS  NOT  ORTHOGONALIZED 


wliere  IW'(n)]  and  [W(n)]  are  M x K for  M regions  or  medium  types,  and 
each  row  of  the  matrices  is  a set  of  weights  to  be  multiplied  over  the  K- 
measurod  event  discriminant  X.  The  detection  measure  D of  M compon- 
ents is  obtained  as: 


D = [ W'(n)  ] X , 


(Iir-18) 


The  resultant  set  of  discriminants  ^ can  be  interpreted  as  projections  on 
unit  vectors,  each  designating  a signal  type.  For  example,  m=  1 for  central 
Asian  events  and  m=  2 for  the  western  United  States  events: 


-1 


1.  1027 
-0. 3364 


-0.  3364 
1,  10  27 


(111-19) 


The  resultant  discriminants  are  plotted  on  Figure  III-5,  These  are  the  orthog- 
onal vector  component  projections  of  the  event  detectability  vectors.  These 
can  be  compared  to  W(n)  • 2£  plotted  on  Figure  III-4.  For  these,  the  vectors 
_U(n)  are  not  orthogonal,  but  are  least  squares  solutions  normalized  to  an 
expected  value  of  1.0  for  presumed  explosions  and  0.0  for  earthquakes.  In 
both  cases,  no  weight  was  given  to  reduction  of  large  negative  deviations. 

The  data  for  Figure  III-4  are  given  under  columns  (5)  and  (6)  of  Table  III-7. 

The  data  for  Figure  III-5  are  given  under  columns  (7)  and  (8)  of  Table  III-7, 

The  data  under  columns  (1)  and  (2)  are  the  discriminants  using  the  average 
detectability  of  presumed  explosions  from  region  1 (central  Asia)  and  from 
region  2 (western  United  States).  Under  columns  (3)  and  (4)  are  the  orthogo- 
nalized  projections.  These  perform  less  effectively  than  the  least  squares 
discriminants  because  of  the  highly  correlated  state  of  the  spectrum  levels 
and  slopes  show'n  on  Table  III-6. 


Also  on  Table  III-7  are  shown  hyper-space  projections  of  the 
discriminants.  Negative  detectabilities  are  weighted  by  zero  and  others  are 
weighted  by  one.  The  or thogonali zed  projections  into  M source  regions  arc 


Ctrntrdl  As\a  Presumed  Explosion  Region 


Western  U.S.  Presumed  Explosion  Region 
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MULTIVARIATE  DISCRIMINANTS  ORTHOGONALI ZED 
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given  as  a length  and  the  square  of  the  direction  cosines.  For  M=  2 (central 
Asia  and  western  United  States)  the  least  squares  discriminant  is  shown  by  a 
length  under  column  (9),  squared  direction  cosine  of  the  region  1 component 
(central  Asia)  under  column  (10),  and  squared  direction  cosine  with  region  2 
under  column  (11). 

The  length  and  squared  direction  cosine  components  ar<>  shown 
for  the  discriminant  derived  from  the  mean  detectabilities  of  region  1 and 
region  2 under  columns  (12),  (13),  and  (14).  For  the  earthquakes  shown  in 
Table  III-7  (page  1),  the  length  of  the  least  squares  projection  was  zero  for 
40  percent  of  the  complex  earthquakes  and  20  percent  of  the  simple  earthquakes. 
For  40  percent  of  each  type  of  earthquake  the  dominant  projection  w'as  toward 
region  1.  For  40  percent  of  the  simple  earthquakes  and  20  percent  of  the  com- 
plex earthquakes  the  projection  was  toward  region  2.  The  lengths  of  the 
earthquake  projections  w'ere  all  less  than  the  length  of  presumed  explosions 
with  the  exception  of  the  tw’O  outliers  shown  on  page  2 of  Table  III-7. 


SECTION  IV 


SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS  | 

1 

The  need  to  perform  interactive  discrimination  processing  in 
the  surveillance  mode  stems  from  the  following  problems: 

• The  methodology  which  optimizes  discrimination  processing  is 

not  yet  known.  Therefore,  programmatic  changes  in  discrim- 
ination processing  procedures  will  need  to  be  rapidly  imple-  | 

mented  on-line.  An  analyst  can  be  computer-a^  sisted  and  i 

prompted  under  the  constraint  of  a programmed  set  of  standard  i 

operating  procedures  for  performing  the  discrimination  process- 
ing. Any  such  programmed  changes  in  the  analyst's  procedures 

must  be  shown  to  significantly  improve  the  analyst's  discrimi- 
nation performance  using  his  old  set  of  standards  as  a baseline. 

• Certain  tasks  such  as  picking  the  start  time  of  an  event,  detec- 
tion and  removal  of  echos,  absorption  correction,  and  other 
interpretive  tasks  may  be  performed  better  by  an  analyst  than 

bv  any  automatic  procedure.  But  the  introductions  of  analyst  ^ 

interpretation,  even  if  constrained  by  standard  operating  pro- 
cedures, trade  off  reproducibility  of  results.  The  projected  ■ 

gains  of  such  processing  must  be  verified  and  shown  to  be  sta- 
tistically significant.  , 

• In  some  difficult  tasks  such  as  discrimination  processing, 
interactive  processing  by  an  analyst  is  a necessary  transition 
to  more  optimal  and  reproducible  automatic  processing.  The 
exact  definition  of  the  task  and  the  criteria  for  optimizing  the 
task  must  be  learned. 


• As  optimum  automatic  processing  procedures  evolve,  the 

analyst  will  function  primarily  in  response  to  emergency  con- 
ditions. Some  of  these  are:  to  detect  and  cope  with  overload- 

ing due  to  event  swarms,  etc.;  to  perform  extended  and  more 
subjective  analyses  on  events  of  doubtful  origin. 

The  development  and  demonstration  of  the  Short-Period  Earth- 
quake/Explosion Discrimination  Package  (SPEED)  partially  fulfills  the  above 
needs  for  interactive  discrimination  processing.  Our  demonstrated  capabili- 
ties have  indicated  that  a substantial  fraction  of  the  required  interactive  pro- 
cessing capability  has  been  developed.  Also,  the  considerable  amount  of  sup- 
port data  processing  for  editing  events,  selecting  channels,  and  for  time 
alignment  of  stations  has  been  developed 

Some  of  the  principal  problems  remaining  to  be  solved  are  as 

follows: 

• Automation  of  support  data  processing: 

select  out  bad  channels 
time  align  channels 
pick  event  start  time 

detect  and  correct  station  or  component  polarity  reversals 

• Extend  and  provide  needed  improvements  to  the;  SPEED  options: 

correct  time  series  for  absorption  with  option  of  repro- 
cessing data  corrected  for  absorption 

provide  analyst  capability'  to  pick  start  times  with  auto- 
matic correction  for  initial  offset 

improve  cepstrum  analysis  to  more  unambiguously  deter- 
mine echo  delays;  if  possible,  automatically  determine 
echo  delays  and  derive  source  parameter  information 
from  the  cepstrum  of  the  signal  waveform 


provide  an  option  to  correct  the  event  waveform  for 
multiple  S3cattered  echos  instead  of  just  for  single  integer 
delay  echos 

change  the  Variable  Frequency  Magnitude  (VFM)  program 
to  measure  niagnitudes  and  rolloffs  at  an  analyst- selected 
set  of  frequencies  instead  of  for  just  a single  high  and 
low  frequency  magnitude  measurement 
expand  the  list  of  short-period  discriminant  variate 
measurements  and  possibly  also  include  long-period  dis- 
criminant variate  measurements  which  improve  overall 
discrimination  performance, 

• Extend  seismic  command  language  capability: 

compile  event  tables  of  discrimination  variates  derived 
from  SPEED  processing 

compile  event  population  plots  of  discrimination  variates 
derived  from  SPEED  processing 

compile  output  tapes  of  standard  format  discriminant 
information  for  extended  least  squares  multivariate  dis- 
criminant analysis. 

These  recommended  improvements  may  not  all  be  achievable 
due  to  storage  limitations  of  the  PDP-15  computer.  In  some  cases  they  may 
have  to  be  performed  as  support  processing  on  the  IIlM  360/44  rather  than  as 
SPEED  options  when  such  storage  limitations  are  eni-ountered. 

In  order  to  establish  a performance  b iseline  on  utilizing  SPEED 
as  a routine  event  discrimination  processor,  a data  ias('  of  representative 
earthquakes  and  explosions  was  processed.  This  baseline  processing  was 
specified  as  simply  as  possible  and  could  be  perforn  ed  reproducibly  in  an 
automatic  mode.  The  baseline  performance  data  was  obtained  by  fixed  pre- 
scribed tapering,  smoothing,  and  inverse  tapering  of  the  ground  disolacement 


measurements.  Several  runs  on  some  of  the  events  indicated  that  the  result- 
ant measurements  were  reproducible.  The  measurements  were  magnitude 
changes  at  specified  high  and  low  frequencies  from  tlie  average  magnitude 
between  1.0  and  1.75  Hz.  It  turned  out  that  the  VFM  method  provides  only 
marginal  discrimination  capability.  It  is  estimated  that  about  half  the  pre- 
sumed explosions  are  detectable  with  a false  alarm  rate  of  about  2-1/2  percent. 
As  would  be  expected,  better  performance  could  be  obtained  by  utilizing  all  of 
the  spectral  values  measured  at  four  low  frequencies  and  four  high  frequencies. 
Unexpectedly,  it  was  found  that  a single  set  of  weights  for  summing  the  dis- 
criminants could  not  detect  presumed  explosions  analyzed  from  central  Asia 
or  from  the  western  United  States.  Two  nearly  orthogonal  sets  of  weights 
were  required.  This  means  that  some  finite  set  of  vectors  of  an  unknown 
dimension  are  necessarily  required  for  a completely  general  multivariate 
spectral  discrimination  of  presumed  explosions.  Even  if  such  a complete 
basis  for  spectral  discrimination  is  derived,  there  is  no  guarantee  that  spec- 
tral discrimination  is  sufficient  to  eliminate  all  false  negatives  associated 
with  certain  types  of  medium  or  geology. 

Th^'re  is  also  no  guarantee  that  certain  regions  may  not  con- 
sistently produce  false  positive  discrimination  of  earthquakes.  The  results 
obtained  for  multivariate  discrimination  indicate  a standard  deviation  of  0,94 
for  the  length  of  positiv'e  deviations.  This  is  a conservative  estimate  ob- 
tained by  neglecting  the  zero  values  of  the  discriminant;  including  them,  the 
standard  deviation  is  0.  86.  The  average  presumed  explosion  discriminant 

. 4 

length  is  3,  21,  Assuming  normal  statistics,  this  would  indicate  a 2.  5 x 1 0 
probability  of  a false  positive  discrimination  at  the  50  percent  level  of  detect- 
ing presumed  explosions.  The  standard  deviation  of  negative  deviations  in  the 
presumed  explosions  population  {possible  false  negatives)  is  0.  53.  The  90  per- 
cent detection  threshold  for  presumed  explosions  is  approximately  2.  5 (aliout 
1.  3 standard  deviations  below  the  mean  of  the  presumed  explosion  population). 


IV-4 


The  probability  of  a false  positive  discrimination  at  the  90  percent  detection 

_3 

level  is  estimated  at  4,  0 x 10  . In  evaluating  the  baseline,  the  two  outliers 

were  obtained,  both  presumed  to  be  peaceful  nuclear  explosions.  From  the 
above  statistics,  these  are  most  likely  not  normal  statistical  variations  but 
indicate  incompleteness  of  the  discrimination  method  in  detecting  every  pos- 
sible kind  of  explosion  event.  Either  additional  vectors  are  needed  for  more 
complete  region  or  source  medium  type  representation  or  more  components 
are  needed  in  the  multivariate  discriminant  vector  X . It  is  probable  that 
more  discriminants  other  than  the  spectral  discriminants  analyzed  are  need- 
ed to  discriminate  all  explosions.  It  is  probable  that  direct  detection  of  focal 
depth  by  methods  more  precise  and  reliable  than  those  presently  used  will  be 
needed  to  eliminate  all  systematic  false  negatives. 
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APPENDIX  A 


EVENT  DATA 

As  part  of  the  quality  control  in  editing  the  NORSAR  data  used 
in  this  report,  a matrix  of  plots  of  all  of  the  event  data  was  prepared  for  use 
by  a data  analyst.  The  plots  included  all  sensors  and  subarray  beams  and 
also  a window  of  noise  preceding  the  event.  The  analyst  by  visual  inspection 
of  this  data  estimated  the  start  time  of  the  event,  eliminated  faulty"  channels, 
channels  with  polarity  reversals,  clipped  channels,  and  those  channels  or 
subarrays  showing  no  visual  indication  of  event  information  on  the  channel 
or  subarray  beam.  The  remaining  channels  were  then  cepstrum  beam- 
formed  and  corrected  for  system  response. 

An  additional  qualiW  control  stop  was  applied  to  the  remaining 
operating  good  channels  by  comparing  the  spectrum  of  the  element  to  a wind- 
ow of  noise  preceding  the  event.  If  at  some  frequency  the  spectrum  of  the 
noise  exceeded  the  spect.-um  of  the  event,  then  that  channel  at  that  frequency 
was  counted  as  containing  no  apparent  event  information.  Eor  each  event,  a 
summary  table  was  prepared  which  at  each  frequency  counted  the  number  of 
channels  with  no  apparent  event  information.  These  are  presented  in  Tables 
A-1  through  A-35  for  all  of  the  events  processed.  The  tables  are  labelled  by 
event  number.  Corresponding  event  description  information  is  shown  on 
Table  III-l.  Plots  of  the  cepstrum  beamform  outputs  are  shown  on  Eigure 
lll-l.  The  column  heading  of  frequency  index  numbers  on  Tables  A-1  through 
A-35  indicate  the  number  of  the  frequency  estimate  of  the  first  row  of  that 
column.  The  index  number  column  on  the  left  side  of  the  table  indicates  the 
frequency  number  increment  starting  from  the  top  of  each  column.  To 
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convert  frequency  index  numbers  to  frequency  in  Hz,  add  the  frequency  index 
number  on  the  left  side  of  the  table  to  the  frequency  index  number  listed  across 
the  top  of  the  table,  subtract  2,  and  multiply  by  0.0098. 

These  tables  give  a powerful  quantitative  check  on  the  amount 
of  useful  event  information  contained  on  all  of  the  channels  used  in  the  cep- 
strum  beamforming  process.  In  case  of  alternative  analyst  picks  of  the 
channels  containing  weak  signals  and  picks  of  the  start  time  of  weak  signals, 
the  tables  containing  lowest  numbers  and  the  broadest  band  of  low  numbers 
indicates  the  preferred  analyst  interpretation.  A quick  method  of  interpreting 
the  tables  is  to  take  the  largest  of  the  first  ten  frequency  numbers  as  the  ex- 
pected number  of  noise  channels.  A positive  indication  of  seismic  information 
is  taken  as  a channel  count  of  noise  at  a given  frequency  number  equal  to  or 
less  than  one-tenth  of  the  largest  count  between  frequency  r.umber  one  and  ten. 
The  number  of  event  indications  by  this  criteria  was  observed  in  the  following 
frequency  number  bands;  1-525,  526-750,  and  751  -1025.  In  the  1 -525  band 
(frequencies  less  than  approximately  5.0  Hz),  all  of  the  events  indicate  sig- 
nificant event  information  using  this  criteria.  In  the  526-750  (frequencies 
approximately  5.  25  Hz  to  7.  5 Hz),  only  one  out  of  20  very  large  earthquakes 
(in  particular,  event  number  22)  indicated  significant  energy  in  this  band  pos- 
sibly due  to  saturation  and  clipping.  Ten  out  of  elev'en  apparent  explosions 
from  eastern  Kazakh  indicated  significant  energy  in  this  band.  The  one  excep- 
tion, event  number  32,  a possible  peaceful  explosion  (PNE),  was  disjilaced 
several  degrees  from  the  other  eastern  Kazakh  epicenters.  In  the  751-1025 
band  (frequencies  greater  than  approximately  7.  5 Hz),  three  out  of  20  c-arth- 
quakes  were  significant  in  this  band;  these  included  ev'ent  number  22  again 
but  also  event  numbers  4 and  87,  both  simple  earthquakes.  Six  out  of  11  of 
the  eastern  Kazakh  events  were  significant  in  t’lls  banti  including  event  num- 
ber 32,  which  was  missofl  in  the  5,  25  Hz  to  7.  5 Hz  barvi. 
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